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Abstract 


The major objective of this study is to develop preliminary concepts for controlling 
orbit, attitude, and structural motions of very large Space Solar Power Satellites (SSPS) 
in geosynchronous orbit. This study focuses on the 1.2-GW “Abacus” SSPS configura- 
tion characterized by a square (3.2 x 3.2 km) solar array platform, a 500-m diameter 
microwave beam transmitting antenna, and an earth-tracking reflector (500 x 700 m). 
For this baseline Abacus SSPS configuration, we derive and analyze a complete set of 
mathematical models, including external disturbances such as solar radiation pressure, 
microwave radiation, gravity-gradient torque, and other orbit perturbation effects. An 
integrated orbit, attitude, and structural control systems architecture, employing electric 
thrusters, is developed. 

A key parameter that characterizes the sensitivity of a satellite to solar radiation 
pressure is the area-to-mass ratio, A/m; the value of A/m for the Abacus satellite is 0.4 
m 2 /kg, which is relatively large when compared to 0.02 m 2 /kg for typical geosynchronous 
communications satellites. Solar radiation pressure causes a cyclic drift in the longitude 
of the Abacus satellite of 2 deg, east and west. Consequently, in addition to standard 
north-south and east- west stationkeeping maneuvers for ±0.1 deg orbit position control, 
active control of the orbit eccentricity using ion thrusters becomes nearly mandatory. 
Furthermore, continuous sun tracking of the Abacus platform requires large control 
torques to counter various disturbance torques. 

The proposed control systems architecture utilizes properly distributed ion thrusters 
to counter, simultaneously, the cyclic pitch gravity-gradient torque, the secular roll 
torque caused by an offset of the center-of-mass and center-of-pressure, the cyclic roll/yaw 
microwave radiation torque, and the solar pressure force whose average value is about 
60 N. A minim u m of 500 ion engines of 1-N thrust level are required for simultane- 
ous attitude and stationkeeping control. When reliability, lifetime, duty cycle, lower 
thrust level, and redundancy of ion engines are considered, this number will increase 
significantly. A significant control-structure interaction problem, possible for such very 
large Abacus platform with the lowest structural mode frequency of 0.002 Hz, is avoided 
simply by designing an attitude control system with very low bandwidth (< orbit fre- 
quency). However, the proposed low-bandwidth attitude control system utilizes a con- 
cept of cyclic-disturbance accommodating control to provide ±5 arcmin pointing of the 
Abacus platform in the presence of large external disturbances and dynamic modeling 
uncertainties. Approximately 85,000 kg of propellant per year is required for simulta- 
neous orbit, attitude, and structural control using 500 1-N electric propulsion thrusters 
with a specific impulse of 5,000 sec. Only 21,000 kg of propellant per year is required 
if electric propulsion thrusters with a specific impulse of 20,000 sec can be developed. 
As I sp is increased, the propellant mass decreases but the electric power requirement 
increases; consequently, the mass of solar arrays and power processing units increases. 
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Chapter 1 

Introduction and Summary 


This chapter provides an executive summary of this report. Detailed technical descrip- 
tions are provided in Chapters 2 and 3. 


1.1 Report Outline 

This report is intended to provide a coherent and unified framework for mathematical 
modeling, analysis, and control of very large Space Solar Power Satellites (SSPS) in 
geosynchronous orbit. 

Chapter 1 presents a summary of major findings and results, and it can be read as 
an executive summary. Chapter 2 provides mathematical models for orbit, attitude, and 
structural dynamics analysis and control design. Chapter 3 presents an integrated orbit, 
attitude, and structural control system architecture, preliminary control analysis and 
design, and computer simulation results. 

Chapters 2 and 3 also contain introductory sections on the basic definitions and 
fundamental concepts essential to mathematical modeling and control of space vehicles. 
These sections summarize many of the useful results in spacecraft dynamics and control, 
and they are primarily based on the material in Space Vehicle Dynamics and Control, 
Wie, B., AIAA Education Series, AIAA, Washington, DC, 1998. 

Chapter 4 provides a summary of the study results and recommendations for future 
research. 

Appendix A begins with a brief description of the general relationship for two-body 
motion, then provides an overview of Encke’s method and how it is carried out in the 
computer program, and ends with a presentation of the expressions used in computing 
the various contributions to the perturbing forces exerted on the two bodies. 
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1.2 Evolution of Space Solar Power Satellites 

A renewed interest in space solar power is spurring a reexamination of the prospects for 
generating large amounts of electricity from large-scale, space-based solar power systems. 
Peter Glaser [1], [2] first proposed the Satellite Solar Power Station (SSPS) concept in 
1968 and received a U.S. patent on a conceptual design for such a satellite in 1973. As 
a result of a series of technical and economic feasibility studies by NASA and DOE in 
the 1970s, an SSPS reference system was developed in the late 1970s. 

The 1979 SSPS reference system, as it is called, featured a very large solar array 
platform (5.3 x 10.7 km) and a gimballed, microwave beam transmitting antenna (1 km 
diameter). The total mass was estimated to be 50 xlO 6 kg. A ground or ocean-based 
rectenna measuring 10 x 13 km would receive the microwave beam on the earth and 
deliver up to 5 GW of electricity. 

In 1995, NASA revisited the Space Solar Power (SSP) concept to assess whether 
SSP-related technologies had advanced enough to alter significantly the outlook on the 
economic and technical feasibility of space solar power. The “Fresh Look” study, [3], 
conducted by NASA during 1995-1997 found that in fact a great deal had changed and 
that multi-megawatt SSP satellites appear viable, with strong space applications. The 
study also found that ambitious research, technology development and validation over 
a period of perhaps 15-20 years are required to enable SSP concepts to be considered 
“ready” for commercial development. 

Recent studies by NASA as part of the SSP Exploratory Research and Technology 
(SERT) program have produced a variety of new configurations of Space Solar Power 
Satellites (SSPS) as reported in Refs. [4] and [5], and shown in Figure 1.1. Some of these 
configurations are based on the passive gravity-gradient stabilization concept. However, 
most other configurations require active three-axis attitude control to maintain contin- 
uous sun tracking of the solar arrays in the presence of external disturbances including 
the gravity-gradient torque. As illustrated in Figure 1.2, a cylindrical configuration, 
which is not affected by the troublesome pitch gravity-gradient torque, has also been 
considered by NASA [6]. 

Two other advanced concepts now under consideration by NASA [6] are shown in Fig- 
ures 1.3 and 1.4. The Integrated Symmetrical Concentrator (ISC) concept, as illustrated 
in Figure 1.3, features ultralight-weight materials and structures which may greatly re- 
duce the projected cost of SSPS. In this concept, mirrors would reflect and focus sunlight 
onto multi-bandgap, thin film photovoltaic arrays located next to a phased-array mi- 
crowave transmitter. On the other hand, the Abacus/Reflector concept, as illustrated 
in Figure 1.4, is characterized by its simple configuration consisting of an inertially 
oriented, 3.2 x 3.2 km solar-array platform, a 500-m diameter microwave beam trans- 
mitting antenna fixed to the platform, and a 500 x 700 m rotating reflector that tracks 
the earth. This study focuses on the 1.2- GW Abacus SSPS configuration. 
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MORPHOLOGY OF SSP CONFIGURATIONS 



Halo 



Clamshell / Sandwich 



(Classifying the Animals in the Zoo) 



Suntower 


Dual Backbone 
Suntower with 
Sub-Arrays 



T - Configurations 


Abacus 



Rigid Monolithic 
Array (1979 
Reference) 


Spin- 
Tensioned 
Monolithic 
Array (Solar 
Disk) 


Figure 1.1: Morphology of various SSPS concepts (Moore [4], [5]). 
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Figure 1.2: A cylindrical SSPS concept with zero pitch gravity-gradient torque (Car- 
rington and Feingold [6]). 
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SERT Systems Integration, Analysis and Modeling 

Integrated Symmetrical Concentrator 


• Massive PMAD estimates motivated 
minimum PMAD configuration 

• Geosynchronous equatorial orbit, mast 
POP 

• Two primary mirror clamshells consisting 
of inflatable flat mirrors (focal length > 10 
km) reflect sunlight onto two centrally- 
located PV arrays, <6X concentration ratio 

• Energy converted to electrical @ PV 
arrays, distributed by cables to 
transmitter, converted to RF and 
transmitted to Earth 

• Clamshells track sun (one rotation per day 
about boom metering structure, seasonal 
beta tilting), transmitter faces rectenna on 
ground 


Concept sized for 1.2 GW 
delivered to Grid on Earth 


3562 x 3644 m 
Clamshells, 

36 456m mirrors 

6378m 

Mast 

Quantum Dot PV 
arrays canted 20° 
846mod ,200mid 

Docking Port 


500 m 
RF Transmitter 



Figure 1.3: Integrated Symmetrical Concentrator concept (Carrington and Feingold [6]). 
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Issues 


Abacus/Reflector Concept 


Prismatic abacus 
frame supports 
lightweight solar 


• In-space construction, assembly, and or deployment of large 
(500m aperture) reflector 

• Surface precision (X/20-X/A0) required by reflector 

• Management of reflector temperature and thermal stresses 

• Azimuth roll-ring and activated links must provide stable platform 
for reflectors 

Benefits 

• Solar collectors always face Sun with very little, if any, shadowing. 

• Solar concentrator uses shifting lens to accommodate seasonal 
beta-tracking, eliminates rotational joints between cells and 
abacus frame. 

• Reflector design eliminates massive rotary joint and slip rings of 
1979 Reference concept. 

• Fixed orbital orientation allows continuous anti-Sun viewing for 
radiators. 

• Abacus structural frame provides runs for PMAD cabling and 
permits “plug and play” solar array approach for assembly and 
maintenance. 

• Triangular truss structure provides reasonable aspect ratio for 
abacus. 

• Activated links provide reflector tilt for target latitude accessability 

• Reduced rotational mass since rotating reflector structure can be 
made much lighter than large planar transmitter array 


concentrators, 

provides structure Arrays of 
for robotic lightweight 

maintenance \ Fresnel 

\ concentrators 


‘k 

iii§ 


"K, 

k. 

issdlfe 


I T' k 


Azimuth roll-ring - 
provides 

E “7 

Ball-screw * 
activated links 
provide reflector tilt 
for various latitude 
pointing 


SMI 


Radiators 
mounted on back 
^ of transmitter 

Orbit 

Normal 


Figure 1.4: 1.2-GW Abacus/Reflector concept (Carrington and Feingold [6]). 
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1.3 Research Objectives and Tasks 

The major objective of this study is to develop preliminary concepts for controlling orbit, 
attitude, and structural motions of very large Space Solar Power Satellites (SSPS) in 
geosynchronous orbit. This study focuses on the 1.2-GW Abacus SSPS configuration 
shown in Figure 1.4. 

The research objectives and tasks of this study, in support of the SSP Exploratory 
Research and Technology (SERT) program of NASA, are as follows: 

• Develop concepts for orbit, attitude, and structural control of very large Space 
Solar Power Satellites (SSPS) using a variety of actuators such as control moment 
gyros, momentum wheels, and electric propulsion thrusters 

• Develop mathematical models, define a top-level control system architecture, and 
perform control systems design and analysis for a baseline Abacus SSPS configu- 
ration in geosynchronous orbit 

• Determine the required number, size, placement, mass, and power for the actua- 
tors to control the orbit, attitude and structural motions of the baseline Abacus 
satellite. Also determine top-level estimates of attitude control system mass and 
propellant consumption per year 

• Further explore advanced control technology toward achieving the mission require- 
ments of future large space vehicles, and provide the space systems designer with 
options and approaches to meet the mission requirements of very large SSPS-type 
platforms 

Because of the limited scope of this study, the following important topics are not 
studied in detail: 

• Thermal distortion and structural vibrations due to solar heating 

• Structural distortion due to gravity-gradient loading 

• Autonomous stationkeeping maneuvers 

• Simultaneous eccentricity and longitude control 

• Attitude control during the solar eclipses 

• Orbit and attitude control during assembly 

• Attitude and orbit determination problem 

• Reflector tracking and pointing control problem 

• Electric propulsion systems 

• Backup chemical propulsion systems 

However, some of these important issues are discussed briefly at appropriate places 
throughout this report. 
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1.4 Abacus SSPS Configuration 

This study focuses on the 1.2-GW Abacus SSPS concept, characterized by a huge (3.2 x 
3.2 km) solar-array platform, a 500-m diameter microwave beam transmitting antenna, 
and a 500 x 700 m rotating reflector that tracks the earth. As illustrated in Figure 1.4, 
some unique features of the Abacus satellite relative to the 1979 SSPS reference system, 
can be described as follows: 

• The transmitting antenna is not gimballed; instead, an azimuth roll-ring mounted, 
rotating reflector provides earth pointing of the microwave beam. 

• The rotating reflector design thus eliminates massive rotary joint and slip rings of 
the 1979 SSPS reference concept. 

• Ball-screw activated links provide reflector tilt for various latitude pointing. 

1.4.1 Geometric Properties 

The three major parts of the Abacus satellite and their dimensions are shown in Figure 
1.5; the mass of each part is given in Table 1.1, together with the total mass and area of 
the spacecraft. The mass of the reflector is approximately 3% of the total mass; therefore, 
the reflector’s mass can be neglected in the analysis of attitude motion, simplifying the 
task in two important respects. First, the Abacus satellite can be treated as a single 
body rather than a multibody spacecraft. When the Abacus satellite is regarded as rigid, 
as is the case in Sec. 2.3, the spacecraft’s moments and products of inertia for a set of 
axes fixed in the solar array do not vary with time. Second, when the unsymmetrical 
mass distribution of the reflector is left out of account, the principal axes of inertia of 
the spacecraft with respect to the spacecraft’s mass center are parallel to the roll, pitch, 
and yaw axes illustrated in Figure 1.5. The moments of inertia for these axes, henceforth 
considered to be principal moments of inertia, are given in Table 1.1. 

The center of pressure is located 100 m below the geometric center of the square 
platform, and the center of mass is located 300 m below the geometric center along the 
pitch axis. The mass of individual components, in units of metric tons, can be found in 
Figure 1.6. 

1.4.2 External Disturbances 

External disturbances acting on the Abacus satellite include: solar radiation pressure 
force, microwave radiation force, gravity-gradient torque, and other orbit perturbation 
forces. Some of these disturbances are summarized in Table 1.2. Disturbance torques 
in units of N-m due to solar pressure, microwave radiation, cm-cp offset, and cm/cp 
location uncertainty can be expressed along the principal axes as: 
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Table 1.1: Geometric and mass properties of the 1.2-GW Abacus satellite 


Solar array mass 

Transmitting antenna mass 

Reflector mass 

Total mass 

Platform area 

Area-to-mass ratio 

Roll inertia 

Pitch inertia 

Yaw inertia 

cm-cp offset 

cm-cp offset (uncertainty) 


21 x 10 6 kg 
3 x 10 6 kg 
0.8 x 10 6 kg 
m = 25 x 10 6 kg 
A = 3200 m x 3200 m 
A/m = 0.4 m 2 /kg 
Ji = 2.8 x 10 13 kg-m 2 
J 2 = 1.8 x 10 13 kg-m 2 
J3 = 4.6 x 10 13 kg-m 2 
200 m (along pitch axis) 
±20 m (along roll axis) 



3200 m 


Figure 1.5: Baseline 1.2-GW Abacus satellite. 
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Abacus/Reflector Concept 

ENTECH Concentrators, Solid State, 1200 MW Delivered from GEO 


System Element 

Mass (MT) 

Comments 

RF Transmitter Array 


Devices, Structure; Input Power = 3614 MWe 

Transmitter Elements 

1156 

Diameter = 500 meters; 83903 Thousand Solid State Devices 

Transmitter Planar Array 

1612 

Mass = 8.21 kg/m2 

Transmitter Array Structure 

281 

Composite Truss Structure @ 1 .43 kg/m2 

Reflector and Bearing Struct 

844 

Only Applicable to Reflector Concepts 

Transmitter Thermal Control 

0 

Integrated, Some TC Included in Transmitter Element Mass 

Add'l Structure Allowance 

117 

Allowance = 3% 

Solar Conversion 


SLA, 1 wing(s) with array dimensions = 1772 m x 6200 m 

Solar Concentrators/Arrays 

4317 

Unit Height = 40 m, Width = 200 m, Mass = 3.664 MT, Power = 3.346 MW 

Add'l Structure Allowance 

129 

Allowance = 3% 

Telecomm & Command 

3 

One set per solar array node (38 sets) 

Add'l Structure Allowance 

0 

Allowance = 3% 

Integrating Structure 

3563 

Abacus, Total Length = 1772 meters 

PMAD 


Cabling & Power Conversion, SPG Power = 3941 MWe; Advanced PMAD 

Cabling 

173 

Total Length = 3162 km @ 0.055 kg/m, Voltage = 100 kV 

Array Converter Mass 

3544 

Mass based on 1178 Converters (1000 V to 100 kV), 3.346 MW Power Out 

Transmitter PMAD Mass 

4362 

Mass Includes Voltage Convertors, Switches, Harness & PMAD Thermal , 3.61 GW 

Rotary Joints, Switches, Etc. 

1 

Thruster Switches Only 

Attitude Control/Pointing 


Sensors, Computers, Control Effectors 

Dry Mass 

452 

Thrusters, CMG's, Sensors etc. at each solar collector 

Propellant 

665 

Delta V = 50 m/s per year for 10 years 

Add'l Structure Allowance 

14 

Allowance = 3% 

SEP Propulsion 

831 

LEO-GEO Transfer Stages 

Dry Mass 

3815 

Thruster Power = 50 kWe (8 Hall Thrusters, Isp = 2000 sec) 

Propellant 

7361 

Krypton 

Add'l Structure Allowance 

114 

Allowance = 3% 

Expendable Solar Arrays 

664 

Thin Film Arrays, 500 W/kg 

Payload Mass 

0 


Satellite Launched Mass (MT) 

33187 

ETO Payloads = 40 MT per launch (830 Launches) 

Satellite Orbited Mass (MT) 

25162 

At 35786 km 


Rectenna Diameter (m) 7450 Harvey Feingold 5/26/00 


Figure 1.6: Mass breakdown of Abacus components (Carrington and Feingold [6]). 
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Table 1.2: Solar pressure and microwave radiation disturbances 


Solar pressure force 
Solar-pressure-induced acceleration 
Solar pressure torque (roll) 

Solar pressure torque (pitch) 
Transmitter /reflector radiation force 
Transmitter /reflector radiation torque 


(4.5E-6)(1.3)(A) = 60 N 
(4.5E-6)(1.3)(A/m) = 2.4 xlO -6 m/s 2 
60 N x 200 m = 12,000 N-m 
60 N x 20 m = 1,200 N-m 
7 N (rotating force) 

7 N xl700 m = 11,900 N-m 


Roll: d\ « 12,000 — 11, 900 cos nt 
Pitch: ^2 ~ 1, 200 
Yaw: d$ ~ — 11, 900 sin nt 

where n is the orbital rate of the Abacus satellite, and t is time. The constant pitch 
disturbance torque of 1,200 N-m is due to the assumed cm-cp offset uncertainty of 20 
m along the roll axis. In addition to these disturbances, gravity-gradient disturbance 
torques are also acting on the Abacus satellite. It is assumed that the electric currents 
circulate in the solar array structure in such a way that magnetic fields cancel out and 
the Abacus satellite is not affected by the magnetic field of the earth. 

1.4.3 Orbit Parameters and Control Requirements 

To describe a satellite orbit about the earth, we often employ six scalars, called the 
six orbital elements. Three of these scalars specify the orientation of the orbit plane 
with respect to the geocentric-equatorial reference frame, also called the Earth-Centered 
Inertial (ECI) reference system, which has its origin at the center of the earth, as shown 
in Figure 1.7. Note that this reference frame is not fixed to the earth and is not rotating 
with it; rather the earth rotates about it. The (A, Y) plane of the geocentric-equatorial 
reference frame is the earth’s equatorial plane, simply called the equator. The X-axis 
is along the earth’s polar axis of rotation. The X-axis is pointing toward the vernal 
equinox , the point in the sky where the sun crosses the equator from south to north on 
the first day of spring. The vernal equinox direction is often denoted by the symbol T. 

The six classical orbital elements consist of five independent quantities which are 
sufficient to completely describe the size, shape, and orientation of an orbit, and one 
quantity required to pinpoint the position of a satellite along the orbit at any particular 
time, as also illustrated in Figure 1.7. 
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Figure 1.7: Orbit orientation with respect to the geocentric-equatorial reference frame, 
also called the Earth-Centered Inertial (ECI) reference system. A near circular orbit is 
shown in this figure. 


Six such classical orbital elements are: 

a = the semimajor axis 
e = the eccentricity 
i = the inclination of the orbit plane 
S2 = the right ascension or longitude of the ascending node 
oo = the argument of the perigee 
M = the mean anomaly 

A traditional set of the six classical orbital elements includes the perigee passage time 
instead of the mean anomaly. The elements a and e determine the size and shape of the 
elliptic orbit, respectively, and t p or M relates position in orbit to time. The angles SI 
and i specify the orientation of the orbit plane with respect to the geocentric-equatorial 
reference frame. The angle oo specifies the orientation of the orbit in its plane. 

Basic orbital characteristics and control requirements for the Abacus satellite in 
geosynchronous orbit are summarized in Table 1.3. 
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Table 1.3: Orbit parameters and control requirements 


Earth’s gravitational parameter 
Geosynchronous orbit (e, i « 0) 
Orbit period 
Orbit rate 
Longitude location 
Stationkeeping accuracy 
Solar array pointing accuracy 
Microwave beam pointing accuracy 


/i = 398,601 km 3 /s 2 
a = 42,164 km 

23 hr 56 min 4 sec = 86,164 sec 
n = 7.292 xlO -5 rad/sec 

TBD 

±0.1 deg (longitude /latitude) 
±0.5 deg for roll/pitch axes 
±5 arcmin 


1.5 Major Findings and Results 

In this section, a summary of major findings and results from this study is presented. 
Detailed technical discussions of the development of mathematical models and a control 
system architecture will be presented in Chapters 2 and 3. 

1.5.1 Technical Issues 

Momentum Storage Requirement 

Assuming that the gravity-gradient torque is the only external disturbance torque acting 
along the pitch axis, we consider the pitch equation of motion of the Abacus satellite in 
geosynchronous orbit given by 

J 2 0 2 = ~y(Js - Ji) sin 20 2 ± u 2 (1) 

where J\, J 2 , and J3 are, respectively, the roll, pitch, and yaw principal moments of 
inertia; 0 2 is the pitch angle measured from the LVLH (local vertical and local horizontal) 
reference frame; n is the orbit rate; and u 2 is the pitch control torque. 

For continuous sun pointing of the Abacus satellite with 0 2 = nt, the pitch control 
torque required to counter the cyclic gravity-gradient torque simply becomes 

372^ 

u 2 = — — (J 3 - Ji)sin2nf (2) 

with peak values of ±143,000 N-m. If angular momentum exchange devices, such as 
momentum wheels (MWs) or control moment gyros (CMGs), are to be employed for 
pitch control, the peak angular momentum to be stored can then be estimated as 

tfmax = y (Js ~ Ji) = 2 x 10 9 N-m-s (3) 
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Table 1.4: A large single-gimbal CMG 


Cost 

$1M 

Momentum 

7,000 N-m-s 

Max torque 

4,000 N-m 

Peak power 

500 W 

Mass 

250 kg 

Momentum / mass 

28 N-m-s/kg 


Table 1.5: A space-constructed, large-diameter momentum wheel [7] 


Momentum = 4 x 10 s N-m-s 
Max torque = 30,000 N-m 
Material = aluminum 
Natural frequency = 0.22 Hz 
Momentum/mass = 66,000 N-m-s/kg 


Rim radius = 350 m 
Mass = 6000 kg 
Max power = 19 kW 
Max speed = 6 rpm 
cost = TBD 


This is is about 100,000 times the angular momentum storage requirement of the In- 
ternational Space Station (ISS). The ISS is to be controlled by four double-gimballed 
CMGs with a total momentum storage capability of about 20,000 N-m-s. The double- 
gimballed CMGs to be employed for the ISS have a momentum density of 17.5 N-m-s/kg, 
and future advanced flywheels may have a larger momentum density of 150 N-m-s/kg. 
Basic characteristics of a large single-gimbal CMG are also summarized in Table 1.4. 

Based on the preceding discussion, it can be concluded that a traditional momentum 
management approach using conventional CMGs (or even employing future advanced 
flywheels) is not a viable option for controlling very large Space Solar Power Satellites. 

To meet the momentum storage requirement of very large SSPS, a concept of con- 
structing large-diameter momentum wheels in space has been studied in the late 1970s 

[7] . An example of such space-assembled, large-diameter wheels is summarized in Table 
1.5. About 5 to 7 such large-diameter momentum wheels are required for the Abacus 
satellite. The concept of lightweight, space-assembled (or deployable, inflatable) large- 
diameter momentum wheels merits further study, but is beyond the scope of the present 
report. 

In an attempt to resolve the angular momentum storage problem of large sun-pointing 
spacecraft, a quasi-inertial sun-pointing, pitch control concept was developed by Elrod 

[8] in 1972, and further investigated by Juang and Wang [9] in 1982. However, such a 
“free-drift” concept is not a viable option for the Abacus satellite because of the large 
pitch attitude peak error of 18.8 deg and its inherent sensitivity with respect to initial 
phasing and other orbit perturbations. 
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Because the pitch gravity-gradient torque becomes naturally zero for cylindrical, 
spherical or beam-like satellites with Ji = J 3 , a cylindrical SSPS configuration was also 
studied by NASA (see Figure 1.2) to simply avoid such a troublesome pitch gravity- 
gradient torque problem. 


Solar Radiation Pressure and Large Area-to-Mass Ratio 

Despite the importance of the cyclic pitch gravity-gradient torque, this study shows 
that the solar radiation pressure force is considerably more detrimental to control of the 
Abacus satellite (and other large SSPS) because of an area-to-mass ratio that is very 
large compared to contemporary, higher-density spacecraft. 

The significant orbit perturbation effect of the solar pressure force on large spacecraft 
with large area-to-mass ratios has been investigated by many researchers in the past [10]- 
[15]. A detailed physical description of the solar radiation pressure can be found in a 
recent book on solar sailing by Mclnnes [14]. The solar pressure effects on formation 
flying of satellites with different area-to-mass ratios were also recently investigated by 
Burns et al. [15]. 

For typical geosynchronous communications satellites, we have 


Area-to-mass ratio A/m ~ 0.02 m 2 /kg 
Solar pressure perturbation acceleration « 0.12 x 10~ 6 m/s 2 
3 tt( 4.5 x 10 -6 )(1.3)A/m 


Ae 


n 2 a 


4.9 x 10 6 per day 
2 


Earth’s gravitational acceleration = 0.224 m/s 
Earth’s oblateness J 2 perturbation = 2.78 x 10~ 6 m/s 2 
Solar gravitational perturbation < 4 x 10~ 6 m/s 2 

2 

Lunar gravitational perturbation < 9 x 10 m/s 
=> Stationkeeping AV « 50 m/sec per year 

AV' 


Am = m 1 — exp 


9l 


sp 


17 kg/year 


where m is assumed as 1000 kg, g = 9.8 m/s 2 , I sp = 300 sec 


For the Abacus satellite, however, we have 

Area-to-mass ratio A/m ~ 0.4 m 2 /kg 
Solar pressure force « 60 N 

Solar pressure perturbation acceleration « 2.4 x 10 -6 m/s 2 
Earth’s gravitational acceleration = 0.224 m/s 2 
Earth’s oblateness J 2 perturbation = 2.78 x 10 -6 m/s 2 
Solar gravitational perturbation < 4 x 10 -6 m/s 2 
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Lunar gravitational perturbation <9x10 6 m/s 


Ae = 


3tt( 4.5 x 10- 6 )(1.3 )A/m 


n 2 a 


1 x 10 per day 


Longitude drift AA = 2Ae « 0.0115 deg/day 
maximum Ae « 0.018 at the mid year 
maximum AA = 2Ae « 2 deg; maximum Aa « 1.8 km 
Orbit eccentricity control is necessary 
(60) (24 x 3600 x 365) 


Am 


5000 x 9.8 


40, 000 kg/year with I sp = 5000 sec 


Typical north-south and east-west stationkeeping maneuvers for the Abacus satellite 
will also require 


Am — ml 1 — exp ( — ) ~ 30, 000 kg/year 

V V 9h P J) 

where m = 25 x 10 6 kg, AV = 50 m/s per year, g = 9.8 m/s 2 , and I sp — 5000 sec. 

The results of 30-day simulations of orbital motion of the Abacus satellite, with the 
effects of the earth’s oblateness and triaxiality, luni-solar perturbations, and 60-N solar 
pressure force included, are shown in Figures 1.8 and 1.9. It is worth noting the extent 
to which eccentricity and inclination are perturbed. 

The initial values used in the simulations correspond to a circular, equatorial orbit 
of radius 42164.169 km; therefore, the initial orbital elements are 

a = 42164.169 km 
e — 0 
i = 0 deg 
— 0 deg 
u) = 0 deg 

The epoch used to calculate the solar and lunar positions, as well as the Earth’s orien- 
tation in inertial space, is March 21, 2000. In order to place the spacecraft at an initial 
terrestrial longitude of 75.07 deg (one of the stable longitudes), a true anomaly 9 of 
253.89 deg is used. 

These elements correspond to an initial position and velocity of 

f = -11698.237/ -40508.869 J + OAkm 
v= 2.954/- 0.853 J + OK km/s 
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Number of orbits 


Figure 1.8: Orbit simulation results with the effects of the earth’s oblateness and triax 
iality, luni-solar perturbations, and 60-N solar radiation pressure force. 



















Table 1.6: Electric propulsion systems for the 1.2- GW Abacus satellite 


Thrust, T 

> 1 N 

Specific impulse, I sp = Tj {frig) 

> 5,000 sec 

Exhaust velocity, V e = I sp g 

> 49 km/s 

Total efficiency, 77 = P 0 /P* 

> 80% 

Power /thrust ratio, P/T 

< 30 kW/N 

Mass/power ratio 

< 5 kg/kW 

Total peak thrust 

200 N 

Total peak power 

6 MW 

Total average thrust 

80 N 

Total average power 

2.5 MW 

Number of 1-N thrusters 

> 500 

Total dry mass 

> 75,000 kg 

Propellant consumption 

85,000 kg/year 


Note: T = mV e , P 0 = |mE e 2 = \TV e , P 0 /T = jV e = ideal power/thrust ratio, Pi/T = 
■^V e , I sp = T/{fng ) = V e /g, V e = I sp g where g = 9.8 m/s 2 , fn is the exhaust mass flow 
rate, Pi is the input power, and P 0 is the output power. 

1.5.2 Control Systems Architecture 

The preceding section illustrates the consequences of solar pressure acting on a spacecraft 
with a large area-to-mass ratio. If left uncontrolled, this can cause a cyclic drift in the 
longitude of the Abacus satellite of 2 deg, east and west. Thus, in addition to standard 
north-south and east- west stationkeeping maneuvers for ±0.1 deg orbit position control, 
active control of the orbit eccentricity using electric thrusters with high specific impulse, 
I sp , becomes mandatory. Furthermore, continuous sun tracking of the Abacus satellite 
requires large control torques to counter various disturbance torques. A control systems 
architecture developed in this study utilizes properly distributed ion thrusters to counter, 
simultaneously, the cyclic pitch gravity-gradient torque and solar radiation pressure. 

Electric Propulsion Systems 

Basic characteristics of electric propulsion systems proposed for the Abacus satellite are 
summarized in Table 1.6. Approximately 85,000 kg of propellant per year is required 
for simultaneous orbit, attitude, and structural control using 500 1-N electric propulsion 
thrusters with I sp = 5,000 sec. The yearly propellant requirement is reduced to 21,000 kg 
if an I sp of 20,000 sec can be achieved. As I sp is increased, the propellant mass decreases 
but the electric power requirement increases; consequently, the mass of solar arrays and 
power processing units increases. Based on 500 1-N thrusters, a mass/power ratio of 
5 kg/kW, and a power/thrust ratio of 30 kW/N, the total dry mass (power processing 
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Magnetic field enhances 
ionization efficiency 


™~~ I 

Magnetic rings B 


1 . Xenon propellant 
injectedj 


Anode 



2. Electrons emitted by hollow 
cathode traverse discharge 
and are collected by anode 


5. Ions are electrostatically 
accelerated through engine 
grid and into space at 30 km/s 


Ion beam 


6. Electrons are injected into 
ion beam for neutralization 


Hollow cathode plasma 
bridge neutralizer 


Figure 1.10: A schematic illustration of the NSTAR 2.3-kW, 30-cm diameter ion thruster 
on Deep Space 1 Spacecraft (92-mN maximum thrust, specific impulse ranging from 
1,900 to 3,200 sec, 25 kW/N, overall efficiency of 45-65%). 


units, thrusters, tanks, feed systems, etc.) of electric propulsion systems proposed for 
the Abacus satellite is estimated as 75,000 kg. 

A schematic illustration of the 2.3-kW, 30-cm diameter ion engine on the Deep Space 
1 spacecraft is given in Figure 1.10, which is formally known as NSTAR, for NASA Solar 
electric propulsion Technology Application Readiness system. The maximum thrust level 
available from the NSTAR ion engine is about 92 mN and throttling down is achieved by 
feeding less electricity and xenon propellant into the propulsion system. Specific impulse 
ranges from 1,900 sec at the minimum throttle level to 3,200 sec. 

In principle, an electric propulsion system employs electrical energy to accelerate 
ionized particles to extremely high velocities, giving a large total impulse for a small 
consumption of propellant. In contrast to standard propulsion, in which the products of 
chemical combustion are expelled from a rocket engine, ion propulsion is accomplished 
by giving a gas, such as xenon (which is like neon or helium, but heavier), an electrical 
charge and electrically accelerating the ionized gas to a speed of about 30 km/s. When 
xenon ions are emitted at such high speed as exhaust from a spacecraft, they push the 
spacecraft in the opposite direction. However, the exhaust gas from an ion thruster 
consists of large numbers of positive and negative ions that form an essentially neutral 
plasma beam extending for large distances in space. It seems that little is known yet 
about the long-term effect of such an extensive plasma on geosynchronous satellites. 
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Orbit, Attitude, and Structural Control Systems 

A control systems architecture developed in this study is shown in Figure 1.11. The 
proposed control systems utilize properly distributed ion thrusters to counter, simul- 
taneously, the cyclic pitch gravity-gradient torque, the secular roll torque caused by 
cm-cp offset and solar pressure, the cyclic roll/yaw microwave radiation torque, and the 
solar pressure force whose average value is 60 N. A control-structure interaction prob- 
lem of the Abacus platform with the lowest structural mode frequency of 0.002 Hz is 
avoided simply by designing an attitude control system with very low bandwidth (< 
orbit frequency). However, the proposed low-bandwidth attitude control system utilizes 
a concept of cyclic-disturbance accommodating control to provide ±5 arcmin pointing 
of the Abacus platform in the presence of large external disturbances and dynamic mod- 
eling uncertainties. High-bandwidth, colocated direct velocity feedback, active dampers 
may need to be properly distributed over the platform. 

Placement of a minimum of 500 1-N electric propulsion thrusters at 12 different 
locations is illustrated in Figure 1.12. In contrast to a typical placement of thrusters 
at the four corners, e.g., employed for the 1979 SSPS reference system, the proposed 
placement shown in Figure 1.12 minimi zes roll/pitch thruster couplings as well as the 
excitation of platform out-of-plane bending modes. A minimum of 500 ion engines of 1- 
N thrust level are required for simultaneous attitude and stationkeeping control. When 
reliability, lifetime, duty cycle, lower thrust level, and redundancy of ion engines are 
considered, this number will increase significantly. 

1.5.3 Attitude and Orbit Control Simulation Results 

Computer simulation results of a case with initial attitude errors of 10 deg in the pres- 
ence of various dynamic modeling uncertainties (e.g., ±20% uncertainties in moments 
and products of inertia, cm location, principal axes, etc.) are shown in Figure 1.13. 
The proposed low-bandwidth attitude control system utilizing the concept of cyclic- 
disturbance accommodating control meets the ±5 arcmin pointing requirement of the 
Abacus platform in the presence of large external disturbances and dynamic modeling 
uncertainties. Proper roll/pitch thruster firings needed for simultaneous eccentricity and 
roll/pitch attitude control can be seen in Figure 1.14. Nearly linear control forces are 
generated by on-off modulation of individual 1-N thrusters, as can be seen in Figure 1.14. 
The total thrusting force from the roll/pitch thrusters #1 through #4 nearly counters 
the 60-N solar pressure force. 

Orbit control simulation results with the effects of the earth’s oblateness and triax- 
iality, luni-solar perturbations, 60-N solar pressure force, and simultaneous orbit and 
attitude control thruster firings are shown in Figures 1.16 and 1.17. In Figure 1.17, 
Fz is the orbit inclination control force and Fx is the solar pressure countering control 
force. It can be seen that the inclination, eccentricity, satellite longitude location, and 
the Z-axis orbit position are all properly maintained. 
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Feedforward Control 


Torque Command Gravity -Gradient T orque 

u 2c= 3n 2 (J 1 -J 3 )(sin2nt)/2 -3n 2 (Jj- J 3 ) (sin 2 0 2 )/2 



Figure 1.11: An integrated orbit, attitude, and structural control system architecture 
employing electric propulsion thrusters. 
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Thrust force direction 



Roll: 1/3 Pitch: 2/4 Yaw: 5/6/7/8 
Orbit Eccentricity, Roll/Pitch Control: 1/3, 2/4 

E/W and Yaw Control: 9/10/11/12 
N/S and Yaw Control : 5/6/7/8 


Figure 1.12: Placement of a minimum of 500 1-N electric propulsion thrusters at 12 
different locations, with 100 thrusters each at locations #2 and #4. (Note: In contrast 
to a typical placement of thrusters at the four corners, e.g., employed for the 1979 
SSPS reference system, the proposed placement of roll/pitch thrusters at locations #1 
through #4 minimizes roll/pitch thruster couplings as well as the excitation of platform 
out-of-plane bending modes.) 
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1.6 Summary and Recommendations for Future Re- 
search 

1.6.1 Summary of Study Results 

The area-to-mass ratio, 0.4 m 2 /kg, is a key indication of the sensitivity of the Abacus 
satellite to solar radiation pressure. Left unopposed, solar radiation pressure can cause a 
cyclic drift in the longitude of the Abacus satellite of 2 deg, east and west. Consequently, 
in addition to standard north-south and east-west stationkeeping maneuvers, active 
control of the orbit eccentricity using high-/ sp electric thrusters becomes mandatory. 

The proposed control system architecture utilizes properly distributed 500 1-N ion 
thrusters to counter, simultaneously, the cyclic pitch gravity-gradient torque, the secular 
roll torque caused by cm-cp offset and solar pressure, the cyclic roll/yaw microwave 
radiation torque, and the solar pressure force of an average value of about 60 N. In 
contrast to a typical placement of thrusters at the four corners, e.g., employed for the 
1979 SSPS reference system, the proposed placement shown in Figure 1.12 minimizes 
roll/pitch thruster couplings as well as the excitation of platform out-of-plane bending 
modes. A control-structure interaction problem of the Abacus platform with the lowest 
structural mode frequency of 0.002 Hz is avoided simply by designing an attitude control 
system with very low bandwidth (< orbit frequency). However, the proposed low- 
bandwidth control system utilizes a concept of cyclic-disturbance accommodating control 
to provide ±5 arcmin pointing of the Abacus platform in the presence of large external 
disturbances and dynamic modeling uncertainties. 

Approximately 85,000 kg of propellant per year is required for simultaneous orbit, 
attitude, and structural control using 500 1-N electric propulsion thrusters with I sp 
= 5,000 sec; yearly propellant consumption is reduced to 21,000 kg if the thrusters 
have an I sp of 20,000 sec. As I sp is increased, the propellant mass decreases but the 
electric power requirement increases; consequently, the mass of solar arrays and power 
processing units increases. The total dry mass (power processing units, thrusters, tanks, 
feed systems, etc.) of electric propulsion systems of the Abacus satellite is estimated as 
75,000 kg, based on 500 1-N thrusters and a mass/power ratio of 5 kg/kW. The peak 
power requirement is estimated as 6 MW based on the peak thrust requirement of 200 
N and a power/thrust ratio of 30 kW/N. 

1.6.2 Recommendations for Future Research 

The baseline control system architecture developed for the Abacus satellite requires 
a minimum of 500 ion engines of 1-N thrust level. The capability of present electric 
thrusters are orders of magnitude below that required for the Abacus satellite. If the 
xenon fueled, 1-kW level, off-the-shelf ion engines available today, are to be employed, 
the number of thrusters would be increased to 15,000. The actual total number of ion 
engines will further increase significantly when we consider the ion engine’s lifetime, relia- 
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Table 1.7: Technology advancement needs for the Abacus SSPS 



Current 

Enabling 

Electric Thrusters 

3 kW, 100 mN 

30 kW, 1 N 


I sp = 3000 sec 

I sp > 5000 sec 


(5,000-10,000 thrusters) 

(500-1,000 thrusters) 

CMGs 

20 N-m-s/kg 

2,000 N-m-s/kg 


5,000 N-m-s/unit 

500,000 N-m-s/unit 


(500,000 CMGs) 

(5,000 CMGs) 

Space- Assembled 


66,000 N-m-s/kg 

Momentum Wheels 


4 xlO 8 N-m-s/unit 

(300-m diameter) 


(5-10 MWs) 


bility, duty cycle, redundancy, etc. Consequently, a 30-kW, 1-N level electric propulsion 
thruster with a specific impulse greater than 5,000 sec needs to be developed for the 
Abacus satellite if excessively large number of thrusters are to be avoided. 

Several high-power electric propulsion systems are currently under development. For 
example, the NASA T-220 10-kW Hall thruster recently completed a 1,000-hr life test. 
This high-power (over 5 kW) Hall thruster provides 500 mN of thrust at a specific 
impulse of 2,450 sec and 59% total efficiency. Dual-mode Hall thrusters, which can 
operate in either high-thrust mode or high-/ sp mode for efficient propellant usage, are 
also being developed. 

The exhaust gas from an electric propulsion system forms an essentially neutral 
plasma beam extending for large distances in space. Because little is known yet about 
the long-term effect of an extensive plasma on geosynchronous satellites with regard to 
communications, solar cell degradation, contamination, etc, the use of lightweight, space- 
assembled large-diameter momentum wheels may also be considered as an option for the 
Abacus satellite; therefore, these devices warrant further study. The electric thrusters, 
CMGs, and momentum wheels are compared in Table 1.7 in terms of their technology 
advancement needs. It is emphasized that both electrical propulsion and momentum 
wheel technologies require significant advancement to support the development of large 
SSPS. 

Despite the huge size and low structural frequencies of the Abacus satellite, the 
control-structure interaction problem appears to be a tractable one because the tight 
pointing control requirement can be met even with a control bandwidth that is much 
lower than the lowest structural frequency. However, further detailed study needs to 
be performed for achieving the required 5-arcmin microwave beam pointing accuracy 
in the presence of transmitter /reflector-coupled structural dynamics, Abacus platform 
thermal distortion and vibrations, hardware constraints, and other short-term impulsive 
disturbances. 
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Although the rotating reflector concept of the Abacus satellite eliminates massive 
rotary joint and slip rings of the 1979 SSPS reference concept, the transmitter fixed to 
the Abacus platform results in unnecessarily tight pointing requirements imposed on the 
platform. Further system-level tradeoffs will be required for the microwave-transmitting 
antenna design, such as whether or not to gimbal it with respect to the platform, use 
mechanical or electronic beam steering, and so forth. 

The following research topics of practical importance in the areas of dynamics and 
control of large flexible space platforms also need further detailed investigation to support 
the development of large SSPS. 

• Thermal distortion and vibration due to solar heating 

• Structural distortion due to gravity-gradient loading 

• Autonomous stationkeeping maneuvers 

• Simultaneous eccentricity and longitude control 

• Attitude control during the solar eclipses 

• Orbit and attitude control during assembly 

• Attitude and orbit determination problem 

• Reflector tracking and pointing control problem 

• Microwave beam pointing analysis and simulation 

• Space-assembled, large-diameter momentum wheels 

• Electric propulsion systems for both orbit transfer and on-orbit control 

• Backup chemical propulsion systems for attitude and orbit control 
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Chapter 2 


Mathematical Models of Large 
Sun-Pointing Spacecraft 

2.1 Introduction to Orbit Dynamics 

This section provides a summary of the basic definitions and fundamentals in orbital 
mechanics. It also provides the necessary background material for a non-Keplerian 
orbit model with various orbit perturbation effects to be discussed later in this chapter. 
Further detailed discussions of orbital mechanics can be found in Ref. [16]. 

2.1.1 Two-Body System 

Consider two particles P\ and P 2 > of masses m 1 and m 2 , whose position vectors from a 
point fixed in an inertial reference frame are given by Ri and R 2 , respectively, as shown 
in Figure 2.1. Applying Newton’s second law and his law of gravity to each particle, we 
write the equations of motion as 


miRi = d 


Gmirri2 


( 1 ) 



Figure 2.1: Two-body problem. 
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m 2 R 2 = 


Gmim 2 


( 2 ) 


where f = R 2 — Pi is the position vector from P x to P 2 , r = |f|, P. t = d 2 Ri/dt 2 is the 
inertial acceleration of Pi, and G = 6.6695 x 10 _11 N-m 2 /kg 2 is the universal gravitational 
constant. 

Eliminating mi from Eq. (1), and m 2 from Eq. (2), and subtracting the first result 
from the second, we obtain 

ff +^ f=0 ( 3 ) 

where r = d 2 f/dt 2 is the inertial acceleration of P 2 with respect to Pi, r = |f|, and 
p = G(nii + m 2 ) is called the gravitational parameter of the two-body system under 
consideration. Equation (3) describes the motion of P 2 relative to Pi in an inertial 
reference frame and it is the fundamental equation in the two-body problem. 

In most practical cases of interest in orbital mechanics, the mass of the primary 
body is much greater than that of the secondary body (i.e., mi 3> m 2 ), which results 
in //, « Grrii. For example, for a sun-planet system, we have p « Pq = GMq, where 
Pq denotes the gravitational parameter of the sun and M© denotes the mass of the sun. 
Likewise, for an earth-satellite system, we have p « p® = GM , ®, where p$ denotes 
the gravitational parameter of the earth and M e denotes the mass of the earth. It is 
worth emphasizing that, in the two-body problem, the primary body is not inertially 
fixed. The two-body problem must be distinguished from a so-called restricted two-body 
problem in which the primary body of mass m\ is assumed to be inertially fixed. Such 
a restricted two-body problem is often described by central force motion of a particle of 
mass m 2 around the inertially-fixed primary body of mass m\. 


Energy Equation 

The energy equation of the two-body system is given by 

^ _ R 

2 r 


(4) 


where v = |u| = |r|, the constant 8 is called the total mechanical energy per unit mass 
or the specific mechanical energy, v 2 /2 is the kinetic energy per unit mass, and — p/r 
is a potential energy per unit mass. This equation represents the law of conservation of 
energy for the two-body system. 


Angular Momentum Equation 

Defining the angular momentum per unit mass or the specific angular momentum as 

h=fxf=rxv (5) 
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we obtain 

dh j* 

— = 0 or h = constant vector (6) 

dt v ' 

Thus we have the law of conservation of angular momentum for the two-body system. 

Since h is the vector cross product of f and v, it is always perpendicular to the plane 

containing r and v. Furthermore, since h is a constant vector, f and v always remain in 

the same plane, called an orbital plane. Therefore, we conclude that the orbital plane is 

fixed in an inertial reference frame, and the angular momentum vector h is perpendicular 

to the orbital plane. 


Eccentricity Vector 


Taking the post-cross product of Eq. (3) with h , finding an expression for a vector 
whose inertial derivative is equal to the preceding cross product, and then integrating, 
we obtain 


^ r l 1 - 

r x h r = constant vector = pe 

r 


(7) 


where a constant vector e is called the eccentricity vector. Note that the constant vector 
pe can also be written as 


pe = v x h — — r = v x (r x v) — — r 

= ( v 2 — — )r — (r • v)v 


Taking the dot product of Eq. (7) with r, we find 


h 2 — pr - pre cos 6 (8) 

where h = \h\, e = |e|, and 6 is the angle between f and e. The angle 9 is called the true 
anomaly and e is called the eccentricity of the orbit. 


Kepler’s First Law 

Equation (8) can be further transformed into the orbit equation of the form: 


h 2 / p 

r = 

1 + e cos 9 

(9) 

which can be rewritten as 

V 

r = 

1 + e cos 9 

(10) 

where p, called the parameter, is defined as 


h 2 

P = — 
P 

(11) 
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Equation (10) is the equation of a conic section , written in terms of polar coordinates 
r and 6 with the origin located at a focus, with 6 measured from the point on the conic 
nearest the focus. Kepler’s first law states that the orbit of each planet around the sun 
is an ellipse, with the sun at one focus. Since an ellipse is one type of conic section, 
Kepler’s first law follows from this equation. The size and shape of the orbit depends 
on the parameter p and the eccentricity e, respectively. 


Kepler’s Second and Third Laws 


The orbital area, A A, swept out by the radius vector r as it moves through a small angle 
A 6 in a time interval At, is given as 


AA 


r(r AO) 


Then the areal velocity of the orbit, denoted by dA/dt, can be shown to be constant, as 
follows: 

dA A A .. 1 o A 6 1 2 - 1, 

— — = hm — — = hm -r — — = -r 0 - -h = constant (12) 

dt At— >o At At— >o 2 At 2 2 v J 

which is a statement of Kepler’s second law: the radius vector from the sun to a planet 

sweeps out equal areas in equal time intervals. 

The period of an elliptical can be found by dividing the total orbital area by the 

areal velocity, as follows: 


P 


A 


nab 


na 2 \/l — e 2 


= 2n\ — 


dA/dt h/2 ^ua(l-e 2 )/2 V H 


(13) 


where a is the semimajor axis and b is the semiminor axis of an ellipse. This can be 

rewritten as 

A-yp 

P 2 = a 3 

which is, in fact, a statement of Kepler’s third law: the square of the orbital period of 
a planet is proportional to the cube of the semimajor axis of the ellipse. Note that the 
ratio P 2 ja 3 is not constant for all planets because /i = G(M© + m 2 ), where M© is the 
mass of the sun and ra 2 is the mass of the planet. Therefore, the ratio differs slightly 
for each planet. 


Kepler’s Time Equation 

Now we introduce a geometrical parameter known as the eccentric anomaly to find the 
position in an orbit as a function of time or vice versa. 

Consider an auxiliary circle, which was first introduced by Kepler, as shown in Figure 
2.2. From this figure, we have 


a cos E + r cos(n — 6) = ae 


(14) 
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Figure 2.2: The eccentric anomaly E of an elliptic orbit. 


where E is the eccentric anomaly and 9 is the true anomaly. Using the orbit equation 

p a( 1 — e 2 ) 


we rewrite Eq. (14) as 


lTecos# l + ecos0 
e T cos 9 

cos E = 


(15) 


(16) 


1 + e cos 9 

Using the fact that all lines parallel to the minor axis of an ellipse have a foreshort 
ening factor of 6/a with respect to a circle with a radius of a, we obtain 


rsin# = -(asinE) = a\/ 1 — e 2 sin E 

a 

Combining this with the orbit equation, we obtain 

\/l — e 2 sin# 


sinE = 


1 + e cos 9 


Furthermore, we have 


E sin E 

t an — = 

2 IT cos E 


1 — e 9 

tan - 

1 Te 2 


(17) 


(18) 


( 19 ) 


from which E or 9 can be determined without quadrant ambiguity. 
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( 20 ) 


Equation (14) can be rewritten as 

r cos 6 = a( cos E — e) 

Thus, squaring Eqs. (17) and (20) and adding them, we obtain 

r = a(l — ecosE) (21) 


which is the orbit equation in terms of the eccentric anomaly E and its geometrical 
constants a and e. 

The area swept out by the position vector r is 

= ( 22 ) 


where t p is the perigee passage time, (t — t p ) is the elapsed time since perigee passage, 
and A is the constant areal velocity given by Kepler’s third law: 

irab irab ab 


A = 


P 


2ir^JaFJp ^ 


(23) 


This area of the ellipse is the same as the area of the auxiliary circle swept out by the 
vector R, multiplied by the factor bja. Thus we have 

ab 

~2 

(E — e sin E) (24) 


/ u . , b . 1 , _ ae . 

l- i t-t r ) = -^aE- J asmE) 


which becomes 


2 


(f — tp) = E — e sin E 


(25) 


where E is in units of radians. 

Defining the mean anomaly M and the orbital mean motion n, as follows: 


M = n(t — t p ) 



(26) 

(27) 


we obtain 

M = E — e sin E (28) 

which is known as Kepler’s time equation for relating time to position in orbit. 

The time required to travel between any two points in an elliptical orbit can be 
simply computed by first determining the eccentric anomaly E corresponding to a given 
true anomaly 0 and then by using Kepler’s time equation. 

Kepler’s time equation (28) does not provide time values, (t — t p ), greater than one- 
half of the orbit period, but it gives the elapsed time since perigee passage in the shortest 
direction. Thus, for 9 > ir, the result obtained from Eq. (28) must be subtracted from 
the orbit period to obtain the correct time since perigee passage. 
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Vernal Equinox 


Figure 2.3: Orbit orientation with respect to the geocentric-equatorial reference frame, 
also called the Earth-Centered Inertial (ECI) reference system. A near circular orbit is 
shown in this figure. 

2.1.2 Orbital Elements 

In general, the two-body system characterized by Eq. (3) has three degrees of freedom, 
and the orbit is uniquely determined if six initial conditions are specified, three of which 
are associated with f at some initial time, and three of which are associated with v = f. 
In orbital mechanics, the constants of integration or integrals of the motion are also 
referred to as orbital elements and such initial conditions can be considered as six possible 
orbital elements. 

To describe a satellite orbit about the earth, we often employ six other scalars, 
called the six orbital elements. Three of these scalars specify the orientation of the 
orbit plane with respect to the geocentric-equatorial reference frame, often called the 
Earth- Centered Inertial (ECI) reference system, which has its origin at the center of the 
earth. The fundamental plane in the ECI system, which is the earth’s equatorial plane, 
has an inclination of approximately 23.45 deg with respect to the plane of the earth’s 
heliocentric orbit, known as the ecliptic plane. A set of orthogonal unit vectors {/, J, K} 
is selected as basis vectors of the ECI reference frame with (A, Y, Z) coordinates, as 
shown in Figure 2.3. 

Note that the ECI reference frame is not fixed to the earth and is not rotating with it; 
rather the earth rotates about it. The ( X , Y ) plane of the geocentric-equatorial reference 
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frame is the earth’s equatorial plane, simply called the equator. The Z - axis is along the 
earth’s polar axis of rotation. The X-axis is pointing toward the vernal equinox, the 
point in the sky where the sun crosses the equator from south to north on the first day 
of spring. The vernal equinox direction is often denoted by the symbol T. 

The six orbital elements consist of five independent quantities, which are sufficient 
to completely describe the size, shape, and orientation of an orbit, and one quantity 
required to pinpoint the position of a satellite along the orbit at any particular time. 
The six classical orbital elements axe: 

a = the semimajor axis 
e = the eccentricity 
i = the inclination of the orbit plane 
Cl = the right ascension of the ascending node 
oo = the argument of the perigee 
M = the mean anomaly 

A traditional set of the six classical orbital elements includes the perigee passage time, 
t p , instead of the mean anomaly, M. 

The elements a and e determine the size and shape of the elliptic orbit, respectively, 
and t p or M relates position in orbit to time. The angles Cl and i specify the orientation 
of the orbit plane with respect to the geocentric-equatorial reference frame. The angle 
oo specifies the orientation of the orbit in its plane. Orbits with i < 90 deg are called 
prograde orbits, while orbits with i > 90 deg are called retrograde orbits. The term 
prograde means the easterly direction in which the sun, earth, and most of the planets 
and their moons rotate on their axes. The term retrograde means westerly direction, 
which is simply the opposite of prograde. An orbit whose inclination is near 90 deg is 
called a polar orbit. An equatorial orbit has zero inclination. 

The line of nodes does not exist for equatorial orbits with zero inclination and also 
the line of apsides does not exist for circular orbits with zero eccentricity. Because a set 
of orbit equations with such classical orbital elements has a singularity problem when 
e = 0 or sin i = 0, we often employ the so-called equinoctial orbital elements defined in 
terms of the classical orbital elements, as follows [16]: 


a = a 

Pi= e sin(Q + oo) 
P 2 = ecos(f2 + o;) 
Q i = tan(i/2)sinfl 
Q 2 = tan(i/2)cosfl 

£ = Cl + oo + M 

where i is called the mean longitude. 
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Figure 2.4: Perifocal reference frame. 


2.1.3 Orbital Position and Velocity 

Given the geocentric-equatorial (X,Y,Z) reference frame with basis vectors { I, J, K } 
and a perifocal ( x,y,z ) reference frame with basis vectors {i,j, k}, the position vector 
is represented as 

f=XI + YJ + ZK = xi+yj + zk (29) 

The position vector f can also be expressed as 

f = X'T + Y'.f + Z'K' = X"I" + Y"J" + Z"K" (30) 


where {X ' , Y' , Z') and (X", Y", Z") are the components of the position vector r in two 
intermediate reference frames with basis vectors {I',J',K r } and {/", J", K"}, respec- 
tively. 

The perifocal reference frame is then related to the geocentric-equatorial reference 
frame through three successive rotations as follows: 


' r ' 


cos Q sin Q, 0 


’ / ' 

j' 

= 

— sin Q cos f 1 0 


J 

K' 


0 0 1 _ 


K 



I" ' 


' 1 0 o' 



V ' 


J" 

= 

0 cos i sin* 



J' 

K" 


_ 0 —sin* cos* _ 



K' 


r n 





r 


i 


cos a; sin a; 0 



I" 


j 

= 

— sina; cosa; 0 



J" 


k 


0 0 1 



K" 


(31a) 

(31b) 


(31c) 
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The orbital elements Q, i, and u are in fact the Euler angles of the so-called C 3 (w) 
Ci(i) C 3 (fl) rotational sequence. 

By combining the sequence of rotations above, we obtain 



cos a; 
— sin a; 
0 


which becomes 


where 


sin a; 0 
cos a; 0 

0 1 



1 

0 

0 

0 

cosz 

sinz 

0 

— sinz 

cosz 


C n 

C 12 

c lz 

C 21 

C’22 

G 23 

C 31 

C 32 

C 33 


cosfl 
— sin fl 
0 


"I 

‘ I " 


J 

_ 

K 


sin Q 
cosQ 

0 


~ 

’ I ' 


J 

_ 

K 


( 32 ) 


C\\ - cos ii cos u) — sin sin ui cos i 
C\2 = sin Q cos u + cos Q sin u> cos i 
C*i3 = sin uj sin i 

C21 = — cos Q sin uj — sin il cos uj cos % 

C22 = — sin Q sin o' + cos H cos u cos i 

G 23 = cos uj sin i 

C31 = sin Q sin i 

C32 - — cos Q sin i 

C33 = cosi 


The matrix C = [Cy] is called the direction cosine matrix which describes the orientation 
of the perifocal reference frame with respect to the geocentric-equatorial reference frame. 

The components ( x,y,z ) of the position vector in the perifocal reference frame 
are then related to the components ( X , Y, Z ) of the position vector in the geocentric- 
equatorial reference frame via the same direction cosine matrix C as: 


X " 


\C U Cv2 c 13 1 


' X ' 

y 

= 

CO 

CD 

CM 

cd 

CD 


Y 

z _ 


/~r s~i 

_ ^31 ^32 <^33 


Z 


( 33 ) 


Since the direction cosine matrix C is an orthonormal matrix, i.e., C 1 
have 


X ' 


'C n C 21 C 3 1‘ 


X 

Y 

= 

C\2 C 22 C Z 2 


y 

Z 


1 

CO 

CD 

CO 

CD 

CO 

CD 

1 


_ z _ 


The components of the velocity vector represented as 


C T , we also 


( 34 ) 


v = XI + YJ + ZK = xi + yj + zk 


( 35 ) 


41 



are also related as 


X 


C n 

C 12 

Cis ‘ 


y 

= 

C 21 

C 22 

C 23 


z 


Csi 

C 32 

C 33 



X ■ 


C n 

C 21 

C 3 i ‘ 


Y 

= 

C 12 

C 22 

C 32 


Z . 


c 13 

C 23 

C 33 



Using the following relationships 


X 


r cos 6 


" X " 


— y p/psinO 

y 

= 

rsinO 

and 

y 

= 

J Pl lp(e + cos 0) 

V 

z 


0 


z 


o 


where p = a(l — e 2 ) is the parameter, we also obtain 


X ' 


' C n C 2 i Csi ' 


r cos 0 


Y 

= 

Cl2 C22 C 32 


r sin^ 


Z 


_ Cis C23 ^33 


0 


X ' 


'C n C 21 C 31 - 


— \/ h /P sin 0 

Y 

= 

C\2 C 22 C32 


\Jpjp{e-\- cos0) 

Z _ 


_ C13 C23 C33 


0 


( 36 ) 


(37) 


(38) 


(39) 

(40) 


2.1.4 Geosynchronous Orbits 

If the period of a satellite in a circular prograde equatorial orbit is exactly one siderial 
day (23 h 56 min 4 s), it will appear to hover motionlessly over a point on the equator. 
Such a satellite, located at 42,164 km (« 6.6 R®) from the earth center (or at an altitude 
of 35,786 km), is called a geostationary satellite. A satellite with the orbital period of 
one siderial day but with a non-zero inclination is called a geosynchronous satellite. Its 
ground track is often characterized by a “figure-eight” curve. Note that regardless of 
the satellite’s orbital inclination, geosynchronous satellites still take 23 hr 56 min 4 s to 
make one complete revolution around the earth. 


2.2 Orbital Perturbations 

Thus far in this chapter, we have considered two bodies whose relative motion is de- 
scribed by an ideal or Keplerian orbit in which the plane of the orbit is fixed in inertial 
space. The Keplerian orbit is a consequence of the assumptions that the primary body 
has a spherically symmetric mass distribution, the second body is a particle, and the 
only forces exerted on the two bodies are those of mutual gravitational attraction. In 
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general, however, the mass of the primary body is distributed aspherically; the two 
bodies are subject to the gravitational attraction of other bodies, and to other pertur- 
bational forces. As a result, the orbit of the two bodies is non-Keplerian, and the plane 
of the orbit does not remain fixed in inertial space. The small deviations from the ideal 
Keplerian orbital motion are often called orbital perturbations. This section presents a 
non-Keplerian orbit model of satellites influenced by the earth’s oblateness and triaxi- 
ality, gravitational perturbations from the Sun and Moon, and solar radiation pressure 
force. 


2.2.1 Non-Keplerian Orbit Dynamics 

Consider the general equation of motion of a satellite about the earth described by 

f+^f=/ (41) 


where f is the position vector of the satellite from the center of the earth, f indicates the 
second derivative of r with respect to time in an inertial reference frame, p « p®, and 
/, called the perturbing acceleration, represents the resultant perturbing force per unit 
mass acting on the satellite, added to the negative of the resultant perturbing force per 
unit mass acting on the earth. The position of a satellite acted upon by the perturbing 
acceleration is often referred to a plane containing f and r , called the osculating orbital 
plane. 

Taking the dot product of Eq. (41) with f yields 


u . H ^ u 

# I # 


which is rewritten as 

d ( v 2 p\ 

dt \ 2 r) 

Substituting the specific energy £ defined as 



V 2 fJ, p 

2 r 2 a 


into Eq. (43), we obtain 


a = 


2a? 


Note that £ is not a constant unless / = 0 or / • r = 0. 
Taking the cross product of Eq. (41) with r, we have 


f x r = r x f 


(42) 

(43) 


(44) 


(45) 
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( 46 ) 


Differentiating the specific angular momentum defined as 

h = fxf 


we obtain 


h = r x r = fxf 


Note that h is not a constant vector unless / = 0orf x/=0. 
Taking the post-cross product of Eq. (41) with h, we have 


which is rewritten as 


rxh + x h = f x h 


^ (r x h — — f) = r x h + f x h 


Substituting the eccentricity vector e defined as 

^ u r 
/xe = r x h r 

r 

and Eq. (47) into Eq. (49), we obtain 

He = f x (r x /) + / x h 


(47) 

(48) 

(49) 

(50) 

(51) 


Here, e is not a constant vector unless the right-hand side of Eq. (51) is zero. 

Let e r , ee, and e z be unit vectors along the radial vector direction, the transverse 
orbit direction, and the direction normal to the orbit plane, respectively, such that 
e r x e$ = e z . Then the perturbing acceleration / and the velocity vector v = f are 
represented in terms of the unit vectors {e r , e$, e z }, as follows: 


/ — fr&r + fe&e + fz&z (52) 

v = v r e r + v e e e + v z £ z (53) 


We also have 


v r = r = 


esin^ 


Ve = r0 = (1 + ecos0) 


(54) 

(55) 


and v z = 0 due to the assumptions of the osculating orbit. Consequently, the term / • r 
in Eq. (44) becomes 

f ■ r = f r v r + f e vg (56) 


and we obtain 


a = 


2a 2 


{/ r esin# + f s ( 1 + ecos^)} 


(57) 
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Differentiating the specific angular momentum vector expressed as 

h = k 

where k (= e z ) is a unit vector normal to the orbit plane, we obtain 

2 y p 

Furthermore, we have 

k = {QK + if + uk) x k 
= Q sin if" — iJ" 


where I" is a unit vector toward the ascending node and J" is orthogonal to I" (see 
Figure 2.3). Thus, we have 


1 fa j. 


pk + ^p/p(Q sin il” — iJ" 


The term r x / in Eq. (47) is also written as 

f x / = rf$k — rf z e$ 

In terms of unit vectors I", J", and k, this equation becomes 

fx/ = rfek — rf z [— sin(o; + 0)1” + cos {ui + 0) J") 

Since h = r x /, equating the coefficients of Eqs. (61) and (63) gives 


. 'f f 

Q sin i = ■ z sin(o; + 0) 

s/W 

% = J-JjL- cos(u) + 0) 

\[Wp 

Differentiating the relation, p = a( 1 — e 2 ), gives 


e= 2^[a(!-e 2 )-p] (67) 

Combining this equation with Eqs. (57) and (64) and using the following relationships 

p = a(l — e 2 ) 
r = a(l — ecos E) 
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we obtain 


[f r sin 0 + f e ( cos 0 + cos E )] 


Similar to the preceding derivation of e, the rotation of the major axis can also be 
obtained as 

lo + Cl cost = J ^ [— f r cos6 + fg( 1 + r/p) sin 0] (69) 

Since the mean anomaly is defined as M = E — e sin E, we obtain 

M = E — esinE — eEcosE (70) 

and E may be found by differentiating r = a(l — ecosE), as follows: 

• r — a(l — ecosE) + aecosE . > 

E = : ^ (71) 

ae sin E 

where r = Jp/p e sin 0. Combining these relationships and using 


cos E = 


sin E = 


e + cos 0 
1 + e cos $ 
\/l — e 2 sin# 
1 + e cos 0 


we obtain 


♦ 9 T f 1 — 6^ 

M = n j H [f r cos 0 — f$(l + r/p) sin 0] 


na * nae 


In summary, we have the so-called Gauss form of Lagrange’s planetary equations as 


a = -^ = [/ r esin# + /0(l + ecos0)l 

\[W 

e = fr sin# + f e (cos# + e + C0S M 
y p [ \ 1 + e cos 0 ) 

■. _ rfz cos(o ; + 0) 

s/W 

a = rfz sin(o; + 0) 

y/JIp sin i 

f z r sin( lu + 0 ) cos i 1 [p r . . , ... 

w = -- )— ■ . */- [fr cos 0 — f e (l + r/p) sm( 


^/jupsin* 

2 rf r ,1 - e 2 


M=n-^f + 


na * nae 


[fr cos 0 - /<?( 1 + r/p) sin 0] 


where p = a(l — e 2 ), n = Jpfa?, and r = p/( 1 + ecos0). 
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Because this set of equations has a singularity problem when e = 0 and/or sun = 0, 
another set of equations that are free of singularities is often considered by employing a 
set of the so-called equinoctial orbital elements defined in terms of the classical orbital 
elements, as follows [16]: 

a = a 

Pi = e sin(G -|- uj'j 
P 2 = ecos(Q + u) 

Qi = tan(i/2)sinf2 
Q 2 = tan(i/2) cos 
£ = £! + ui + M 

where l is called the mean longitude. 

Furthermore, using the true longitude and eccentric longitude defined, respectively, 
as 


L — £1 + u) + 0 

(79) 

K = Q + lo + E 

(80) 

we rewrite Kepler’s orbit equation, r = a( 1 — cos E), as 


r = o(l — Pi sin K — P 2 cos K) 

(81) 

and Kepler’s time equation, M = E — e sin E, as 


£ = K + Pi cos K — P 2 sin K 

(82) 


The true longitude, L, can be obtained from the eccentric longitude, K, using the fol- 
lowing relationships [16]: 

sin L = - \( 1 ^-rPo) sin K + ~^-rPiP 2 cos K - Pi] (83) 

rLV a + b / a + b 

cos L = - \ ( 1 a -P?) cos K + — t -P 1 P 2 sin K - P 2 ] (84) 

rLV a + b V a + b J vy 

where 

a _ 1 _ 1 

a + b 1 + y/l - e 2 ^1 _ pa _ P 2 

From Battin [16], we have Gauss’ variational equations in terms of the equinoctial ele- 
ments, as 

r t ) 

a = -r- (P 2 sin L — Pi cos L)f r + -f e (85) 

hi r J 
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A 

A 

Qi 

Q2 

i 


r 


fe - P 2 {Qi cos L-Q 2 sin L)f z J 
— sin L / r + P 2 4- ^14 — ^ cos T fo Pi{Qi cos T — Q 2 sin _L)/ Z |- 


h\r COsLfr + 
r fp 

h 

— (1 + Qi 4- Q 2 ) sinL f z 

— (1 + Qi 4- Q 2 ) cosL f z 


1 4- - ) sin L 

r 


r 

n ~h 


(-) (P\ sin L + P 2 cos L) + — 
\rj a 


+ ~j / (l 4- (Pi cos L - P 2 sin L)f e + (Q 1 cos L - Q 2 sin L)/* j 


( 86 ) 

(87) 

( 88 ) 

(89) 

(90) 

(91) 


where 


& = <V 1 - P? - P 2 2 
h = nab 

— = 1 4- Pi sin L 4- P 2 cos L 

r 

r h 

h n(l + Pi sinL 4- P 2 cosL) 

By defining 

f*= X/4- FJ4- ZA 
/= P x / + F y J + F z K 


where ( X , Y,Z) are the so-called ECI coordinates, we also obtain the orbital dynamic 
equations of the following simple form: 


X = 

X „ 

~ ^ + F X 

(92a) 

Y = 

Y 

— 1^—^ + Fy 

(92b) 

Z = 

Z 

~ ^ 

(92c) 


where r = sfX 2 + Y 2 + Z 2 . 


2.2.2 Asphericity 

The earth is not a perfect sphere; it more closely resembles an oblate spheroid, a body 
of revolution flattened at the poles. At a finer level of detail the earth can be thought 
of as pear shaped, but the orbital motion of geosynchronous spacecraft can be analyzed 
adequately by accounting for the mass distribution associated with the polar flattening 
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and ignoring the pear shape. The equatorial bulge caused by the polar flattening is 
only about 21 km; however, this bulge continuously distorts the path of a satellite. 
The attractive force from the bulge shifts the satellite path northward as the satellite 
approaches the equatorial plane from the south. As the satellite leaves the equatorial 
plane, the path is shifted southward. The net result is the ascending node shifts or 
regresses; it moves westward when the satellite’s orbit is prograde, and eastward for 
a retrograde orbit. In this section, we analyze the effects of the earth’s oblateness, 
characterized by the gravitational coefficient J 2 , on the precession of the node line and 
the regression of the apsidal line of satellite’s orbits. 

The equatorial cross-section of the earth is elliptical rather than circular, with a 65 
m deviation from circular; thus, an oblate spheroid (with a circular cross section at the 
equator) is a less precise representation of the earth than an ellipsoid with axes of three 
distinct lengths. When modeling the earth as an ellipsoid, one therefore refers to its 
triaxiality. The tesseral gravitational harmonic coefficient J 22 of the earth is related 
to the ellipticity of the earth’s equator. There are four equilibrium points separated by 
approximately 90 deg along the equator: two stable points and two unstable points. The 
effect of the triaxiality is to cause geosynchronous satellites to oscillate about the nearest 
stable point on the minor axis. These two stable points, at 75° E longitude and 255° E 
longitude, are called gravitational valleys. A geosynchronous satellite at the bottom of 
a gravitational valley is in stable equilibrium. Satellites placed at other longitudes will 
drift with a 5-year period of oscillation; thus, they require “east-west” stationkeeping 
maneuvers to maintain their orbital positions. The stable equilibrium points are used 
among other things as a “junk-yard” for deactivated geosynchronous satellites. 

2.2.3 Earth’s Anisotropic Gravitational Potential 

As discussed in Section 2.2.2, the earth’s shape is better represented by an ellipsoid 
than a sphere, and its mass distribution is not that of a uniform sphere. To account 
for the nonuniform mass distribution and the resulting nonuniformity in the earth’s 
gravitational field, a gravitational potential is given in general terms by a series of 
spherical harmonics, 


,, f °° n / H \ n 'j 

£/ e (r, (/), A) = - < 1 + XT H ( — ) -Pnm (sin 4>)[C nm cos m\ + S nm sinraA] } (93) 

r l n=2m=0 ' r ^ J 


where the point of interest is described by its geocentric distance r, geocentric latitude <j), 
and geographic longitude A measured eastward from the Greenwich meridian, and where 
R© is the mean equatorial radius of the earth, p = /% is the gravitational parameter 
of the earth, C nm and S nm are the tesseral (n ^ m), sectoral (n = m), and zonal 
(m = 0) harmonic coefficients characterizing the earth’s mass distribution, and P nm is 
the associated Legendre function of degree n and order m. 
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The perturbing gravitational potential, U, such that / = (1 + m 2 /mi)Vf/« VU 
becomes the perturbation acceleration, is then defined as 


u = u m -^ 

r 

oo n (R 


CXJ ft / U \ M 

( — ) Pnm (sin 0) [C nm cos m A + S nm sin m A] 

r ^_ 0 w _ n \ T / 


(94) 


n= 2 m=0 

By separating the terms independent of longitude, we find 

oo n / D \ rt \ 

+ J2 J nm ( — ) P n m(sin0)[cosm(A - A nm )] ^ (95) 

n=2m=l \ T / J 

where the zonal harmonic coefficient, J n , is often defined as J n = — C„o (e g-, J 2 = 
+1082.63 x 10~ 6 ), P n is Legendre polynomial of degree n defined as P n = P n 0 , and 


Jnm — V Cnm + ^nm 


\ 1 , — 1 / Pnm 

A«m = — tan ' 


Pi (sin 0) = sin0 
P 2 (sin0) = (3sin 2 0 — l)/2 
P 3 (sin0) = (5sin 3 0 — 3sin0)/2 
P 4 (sin 0 ) = (35sin 4 0 — 30 sin 2 0 + 3)/8 
Pn(sin0) = (1 — sin 2 ^) 1 / 2 
P 2 i(sin0) = 3 sin 0(1 — sin 2 0) 1//2 
P 22 (sin0) = 3(1 - sin 2 0) 

P3i(sin0) = 3(1 — sin 2 0)^ 2 (5sin 2 0 — l)/2 
P 32 (sin0) = 15(1 — sin 2 0) sin 0 
P33(sin0) = 15(1 — sin 2 0) 3 ^ 2 


Note that P n and P nm can be determined from the following formulas: 

C.M = 7 ^-,^ - D” 


Pnmi^P) 


2 n n! dx n 
(1 — x 2 ) m / 2 d n+m . 2 _ 

¥n\ dx n + m ^ X ~ ’ 


A set of numerical values for the coefficients and constants in Eq. (93) is known as a 
gravitational model. The Goddard Earth Model (GEM) T1 is reported by Marsh et al. 
[21], with P e = 6,378.137 km, and /. t© = 398,600.436 km 3 /s 2 . The normalized values of 
gravitational coefficients in Ref. [21] have been unnormalized into the following values 
for J r and C rs , and used to calculate other parameters, J rs and A„: 
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J 2 = 1082.63E-6 J 3 = — 2.5326E-6 J A = -1.6162E-6 J 5 = -0.22812E-6 
C 21 =0 5 2 i = 0 C 22 = 1.57432E-6 S 22 = -0.903593E-6 

C31 = 2.192406E-6 S 3 1 = 0.269593E-6 C 32 = 0.30862E-6 S 32 = -0.211914E-6 

C 33 = 0.100537E-6 S 33 = 0.197057E-6 

J 22 = 1.81520E-6 J31 = 2.20892E-6 J 32 = 0.37437E-6 J 33 = 0.22122E-6 

A 22 = -0.26052 A31 = 0.12235 A 32 = -0.30085 A 33 = 0.366343 

2.2.4 Earth’s Oblateness 

The effects of the earth’s oblateness on the precession of the node line and the regression 
of the apsidal line of satellite’s orbits can now be analyzed considering the perturbing 
gravitational potential of the oblate earth given by 

U = ^ |-^^(3sin 2 0- 1) - -^^(5sin 3 <^> — 3 sin 4>) j (96) 

where r, <f>, E®, /1, J 2 , and J 3 have the same meanings and numerical values as given in 
Sec. 2.2.3. 

As illustrated in Figure 2.5, the angle between the equatorial plane and the radius 
from the geocenter is called geocentric latitude, while the angle between the equatorial 
plane and the normal to the surface of the ellipsoid is called geodetic latitude. The 
commonly used geodetic altitude is also illustrated in Figure 2.5. 

North Pole (X,Y , Z ) 


Altitude 



Figure 2.5: A two-dimensional view of the oblate earth. 


Ignoring higher-order terms, we consider the perturbing gravitational potential due 
to J 2 , as follows: 

C/ = i^S(l-3sin V) (97) 
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Since the geocentric latitude 4> is related to the orbital elements as: 


. , Z r sin(o; + 0) sin ? . . ... . 

sm 0 = — = = sm w + 0) sm ? 

r r 


Eq. (97) is rewritten as 


( 98 ) 


(99) 


Noting that f= re r and dz = r sin(o; + 0)di, we can express the perturbing acceleration 

^ * „ TT dU _ 1 dU ^ 1 dU _ 

VL^-— e r + -^e 0 + — — — — rr— (100) 


\ 1 3sin 2 ? sin 2 (a; + 6) 


Taking the partial derivatives of U with respect to r, 0, and ?, and substituting them 
into Eq. (100), we obtain the radial, transverse, and normal components of /, as follows; 


fz = — 


3/?J 2 .R| 

2r 4 

3/iJ 2 J?| 

2r 4 

3 /?«/ 2 i?| 
2r 4 


■{1 — 3sin 2 ? sin 2 (u; + 0)} 
sin 2 ? sin2(u? + 0) 

• sin 2? sin(o; + 0 ) 


( 101 ) 

( 102 ) 

(103) 


Substituting Eq. (103) into Eq. (65), we obtain the precession of the node line as: 


fi = - 


3/?J 2 /? e CQg 
r 3 y/jlp 


(104) 


Integrating this equation over an entire orbit of period P yields 


AQ = - 


3pJ 2 /? 2 


cos? 


. f p sin 2 (a; + 0) 


dt 


(105) 


where Afi denotes the change of Q over an entire orbit, assuming that changes in other 
orbital elements are second-order terms. (Note that the average rate of change of ? over 
the orbital period is zero.) 

Since the angular momentum h = \h\ can be expressed as 


h = « r 2 (Hcos? + u; + 0) 


we have 


Co + 1 


\[W> 


(106) 

(107) 
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in which the second-order term Qcos i is further neglected. This equation is used to 
change the independent variable t into (oj + 6), as follows: 


d{u + 9) = dt 


Thus, Eq. (105) can be rewritten as 


AQ = - 


SJ 2 R% . f 2w sin 2 (a; + 9) 


V 


cos % 


d(u) + 9) 


3J 2 R% . f 2v 1 - cos 2(a; + 6) 


■cos? 


p 


2r 


d(u> + 9) 


Performing the integration after a substitution of r = p/(l + e cos 0) yields 


AQ = - 


3ttJ 2 RI 

p2 


cos i + higher-order terms 


(108) 


Dividing this by the average orbital period, P : 27r/n, where n = \J p. / a 3 is the orbital 
mean motion, we obtain the average rate of change of D, as follows: 


n 


2p 2 


ncosi 


(109) 


Similarly, assuming that the eccentricity and the semimajor axis of the orbit remain 
unperturbed by the oblateness of the earth to a first-order approximation, we can obtain 
the average rate of change of u>, as follows: 


(jj 


3 J 2 R, 


■© 


2p 2 


n 



( 110 ) 


For geostationary satellites, we have AQ « —4.9 deg/year and Au « 9.8 deg/year due 
to the oblateness of the earth. 


2.2.5 Earth’s Triaxiality 


The earth’s ellipsoidal shape, or triaxiality, was discussed in Section 2.2.2, where it is 
noted that the equatorial cross-section is elliptical with one axis 65 m longer than the 
other. The elliptical nature of the equator is characterized by the tesseral harmonic 
coefficients C 22 , S 22 , C 32 , S 32 , etc. The primary tesseral harmonic is denoted by J 22 , 
which combines C 22 and S 22 , as 


J22 — 



The earth’s elliptical equator gives rise to a gravitational acceleration that causes a drift 
in the longitudinal position of geostationary satellites, which is a major perturbation 
that must be dealt with. 
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There are four equilibrium points separated by approximately 90 deg along the equa- 
tor: two stable points and two unstable points. The effect of the triaxiality is to cause 
geosynchronous satellites to oscillate about the nearest stable point on the minor axis. 
These two stable points, at 75° E longitude and 255° E longitude, are called gravitational 
valleys. A geosynchronous satellite at the bottom of a gravitational valley is in stable 
equilibrium. Satellites placed at other longitudes will drift with a 5-year period of oscil- 
lation; thus, they require “east-west” stationkeeping maneuvers to maintain their orbital 
positions. The stable equilibrium points are used among other things as a “junk-yard” 
for deactivated geosynchronous satellites. 

Ignoring higher-order terms, the perturbing gravitational potential due to tesseral 
harmonics C 22 and S 22 is defined as 

U = '^—-p-{C 2 2 cos 2A + S 22 sin 2 A) cos 2 4> 

= yj C 22 + Sl 2 cos 2(A — A 22 ) cos 2 (j) 

= J 22 cos 2(A - A 22 ) cos 2 4> (111) 

where r is the geocentric distance, A is the geographic longitude, 4> is the geocentric 
latitude, C 22 = 1.574321 x 10“ 6 , S 22 = -0.903593 x 10" 6 , J 22 = 1.815204 x 10“®, and 
A 22 is defined as 

1 / s \ 

A 22 = - tan -1 ( — ^ ) = —0.26052 rad = —15 deg 

2 \c 22 y 

Using the following relationships 

X = r cos (p cos A 
Y = rcos^sin A 
Z = r sin 0 

we obtain 

U = [C 22 (A 2 - Y 2 ) + 2 S22XY] (112) 

and 


/ = VU = ^S[-5 {C 22 (X 2 - Y 2 ) + 2 S 22 XY}e r 
+2 r(C 22 X + S 22 Y)e 1 + 2 r(S 22 X - 

= fr&r + fe&e + fzk ( 113 ) 

where mutually perpendicular unit vectors e\ and e 2 are fixed in the earth: e\ lies 
in the equatorial plane parallel to a line intersecting earth’s geometric center and the 
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Greenwich meridian, and e 2 lies in the equatorial plane 90 deg eastward of e). Note that 
in Eqs. (112) and (113), X, Y, and Z mark the position of the satellite in an earth-fixed 
coordinate system; they do not have the same meanings as in Eq. (29). 

The radial, transverse, and normal components of / can then be found as 

f r = _ (C '22 cos 2 A + £22 sin 2A) cos 2 (p 

fe = — (C '22 sin 2A — 1 S 22 cos 2A) cos 4> (114) 

fz — ~ (C '22 cos 2 A + iS '22 sin 2A) sin (p cos (p 

For a geosynchronous satellite with and r = a (i.e., <f> = 0), we have 

fr = - (C '22 cos 2A + S 22 sin 2A) 

fe = - (i ^ e (C '22 sin 2 A - S 22 cos 2 A) (115) 

fz = 0 


Using 

^ = - fe where n = J/i/oP 
at n v 

we can express the longitudinal perturbation acceleration as 


3nda 
2 a dt 


and we find 


A = sin 2 A — S 22 cos 2 A) 

a 5 

= 1 j 22 s i n 2(A _ a 22 ) 

a 5 


(116) 

(117) 

(118) 


where A 22 = —0.26052 rad = —15 deg. The equilibrium longitudes, denoted as A*, for 
A = 0 can be found as: 


tan 2A* = ^ 
C 22 

A* = 


-0.90359 x 10" 6 
“ 1.57432 x 10- 6 = 
75, 165, 255, 345 deg 


-0.57395 


and A 22 corresponds to A* of 345 deg. It can be shown that A* of 75 deg and 255 
deg are stable equilibrium longitudes and that A* of 165 deg and 345 deg are unstable 
equilibrium longitudes. 
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Since the stable equilibrium points are separated by ±90 deg from A 22 , we obtain 
the following: 

A = ^^J 22 sin 2(A- A 22 ) 
a 5 

= 22 sin2(A-(A,±7r/2)) 

= _ l8 ^ e J 22 sin 2(A - A s ) 
a 5 

= —0.0017 sin 2(A — A s ) deg/day 2 (119) 

where A s = 75.3 and 255.3 deg (stable longitudes). 


2.2.6 Luni-Solar Gravitational Perturbations 

The gravitational forces exerted by the Sun and Moon on the two bodies of interest, the 
earth and a geostationary satellite, are referred to as the luni-solar perturbation. The 
equation describing motion of a satellite subject to perturbations is given by 

f+ / ( 120 ) 

where / is the perturbing acceleration caused, in this case, by the luni-solar gravitational 
effects on the satellite and earth, described by 


f 


—p>® 





( 121 ) 


and 


f = position vector of satellite from the earth 
= Q + f® for the earth-moon-satellite system 
= R + f© for the earth-sun-satellite system 
Q = position vector of the moon from the earth, Q = 3.84398 x 10 8 m 

R = position vector of the sun from the earth, R = 1 AU = 1.496 x 10 11 m 

f® = position vector of satellite from the moon 
r© = position vector of satellite from the sun 
/% = the earth’s gravitational parameter = 398, 601 km 3 /s 2 

p® = the moon’s gravitational parameter = 4,902.8 km 3 /s 2 

p Q = the sun’s gravitational parameter = 1.32686 x 10 11 km 3 /s 2 


56 



Defining 


f=XI + YJ + ZK 

Q = Qxl + QyJ + QzK 
R = Rxl RyJ Y RzK 

f = F x I + F Y J + F z K 

where (X, Y, Z ) are the ECI coordinates, we can obtain the components of the luni-solar 
perturbation vector, as follows: 

F x ~^l-X + ^ cos#® + ^-y(5 cos 2 #® - 1) (<2*-X)j 

+ ^|-X + | cos#© + ^(5 cos 2 #© — 1) (fl*-X)j (122) 

Fy«||j-y + ^ cos#© + ^(5 cos 2 #® — 1) (Qy-F)} 

+ /f {~ F + f cos# q + ^(5cos 2 #©-1) (Ry-Y) | (123) 

cos#© + ^(5 cos 2 #® — 1) (Qz-Z) j 

+ cos#© + — (5 cos 2 #© — 1) (i?z^^)| (124) 

where #© is the angle between the earth-satellite line and the earth-moon line, #© is the 
angle between the earth-satellite fine and the earth-sun line, and 

Qx = Q{ cos D© cos w©t — sin fl© cos i© sin a ;®t) 

Qr : Q(sin D© cos a;©t+ cos Q® cos «© sin a;®t) 

Qz = Qsini® sinu;®f 

Rx = R( cos D© cos a>©£ — sin f2© cos i© sin u;©£) 
f?y = f?(sin D© cos u Q t+ cos D© cos £ © sin a;©£) 

= R sin i Q sin u Q t 

ujr/j = orbit rate of the moon = 27t/27.3 rad/day 
H® = the right ascension of the moon 
ir/j = the inclination angle of the moon 
a /© = orbit rate of the sun = 27t/365.25 rad/day 
fl© = the right ascension of the Sun 
i Q = the declination angle of the Sun 
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The luni-solar gravitational perturbations for typical geosynchronous communica- 
tions satellites with i « 0 are summarized by Agrawal [22], as follows: 

Lunar gravitational perturbation < 9 x 10 -6 m/s'* 

Solar gravitational perturbation < 4 x 10" b m/s 2 


dt 


3/%r 2 

4/ir| 


sin(Q — Q@) sin i® cos i® 


3 /.i©r 2 


0.478 to 0.674 deg/year 


4hi% 
0.269 deg/year 


sin il sin i Q cos *© 


where Q is chosen as 90 deg, Q® - 0 for min/max i®, and 

/i® = 4.9028 x 10 3 km 3 /s 2 
Ho = 1.32686 x 10 11 km 3 /s 2 
r® « 3.844 x 10 5 km 
r© ~ 1.49592 x 10 8 km 
i® - 18.3° to 28.6° 
i Q = 23.45° 
r = 42, 164 km 
h 129, 640 km"/s 


2,2.7 Solar Radiation Pressure 

The significant orbital perturbation effect of the solar pressure force on large spacecraft 
with large area-to-mass ratios has been investigated by many researchers in the past, 
[10]-[15]. A detailed physical description of the solar radiation pressure can be found in 
a recent book on solar sailing by Mclnnes [14]. The solar pressure effects on formation 
flying of satellites with different area-to-mass ratios were also recently investigated in 
Ref. [15]. 

The solar radiation forces are due to photons impinging on a surface in space, as 
illustrated in Figure 2.6. It is assumed that a fraction, p s , of the impinging photons is 
specularly reflected, a fraction, pd, is diffusely reflected, and a fraction, p a , is absorbed 
by the surface. And we have 

As I Pd I Pa I (125) 

The solar radiation pressure (SRP) force acting on a flat surface is then expressed as 

F = PA{n ■ s) | (p a + Pd)s + ( 2 p s (n ■ s) + ^p d ) ft J (126) 

w r here P = 4.5 x 10”"' 3 N/m 2 is the nominal solar radiation pressure constant, A is the 
surface area, ft is a unit vector normal to the surface, and s is a unit vector pointing 
from the sun to satellite. 
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Sun 


► 


► 


► 



Figure 2.6: Solar radiation pressure force acting on an ideal flat surface (a case with 
45-deg pitch angle (j> is shown here). 

For an ideal case of a perfect mirror with p,i = p a ■ 0 and p s = 1 , we have F t = 0 
and 

F — F n = 2PAcos 2 (j) n 

Also for an ideal case of a black body with p s = pd - 0 and p a = 1, we have 

F = P(A cos 4>) s 

where A cos <p is called the projected area of the surface under consideration. 

For most practical cases of satellites with small pitch angles, the SRP perturbation 
force per unit mass is simply modeled as 

/ - P(1 + p)(A/m)s (127) 

where p is the overall surface reflectance (0 for a black body and 1 for a mirror) and 
A/m is the area-to-mass ratio. 

Defining / f r e r + feee + .fz&z and ignoring the effects of seasonal variations of the 
sun vector, we have 

fr « / sin 0 
fe « / cos 9 

where / = P(1 + p)(A/m). 

From the orbit perturbation analysis, we have 

= —;=={/,. e sin 0 + fe{ 1 + ecos6>)} 

(If x /\ — g2 

-77 = —{fr sin 9 + f e { cos 9 + cos E)} 


59 



For geosynchronous satellites with e « 0, we obtain 

da 2 2/ . 

-77 = ~fe = — sm( 
at n n 

=>• A a = 0 per day 


(128) 


and 


dp 1 

— = — (f r sin 6 + 2 f e cos 0) 
dt na 

= — (/ sin 2 0 + 2/ cos 2 0) 

na 

= — (- + - cos 20 

na V 2 2 

3 tt/ 

— r— per day 
n z a 


Ae 


(129) 


The solar radiation pressure effect on the longitude change can also be found as 

■■--fe 

a 


■■ dn 3 n da 3n 2 

dt 2a dt 2 an J 


= -2£ w 


a 


(130) 


2.2.8 Orbit Simulation Results 

Orbit simulation results for the Abacus satellite with the effects of the earth’s oblateness 
and triaxiality, luni-solar perturbations, and 60-N solar pressure force are shown in 
Figures 2.7 and 2.8. The significance of the orbital perturbation effects on the eccentricity 
and inclination can be seen in these figures. 

Orbit control simulation results with the effects of earth’s oblateness and triaxiality, 
luni-solar perturbations, 60-N solar pressure force, and simultaneous orbit and attitude 
control thruster firings are shown in Figures 2.9 and 2.10. In Figure 2.9, Fz is the 
orbit inclination control force and F x is the solar pressure countering force resulting 
from countering the pitch gravity-gradient torque. It can be seen that the inclination, 
eccentricity, satellite longitude location, and the Z-axis orbital position are all properly 
maintained. The feasibility of using continuous (non-impulsive) firings of ion thrusters 
for simultaneous eccentricity and inclination control is demonstrated. 

The initial values used in the simulations correspond to a circular, equatorial orbit 
of radius 42164.169 km; therefore, the initial orbital elements are 

a = 42164.169 km 
e = 0 
i = 0 deg 
= 0 deg 
oj = 0 deg 
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The epoch used to calculate the solar and lunar positions, as well as the Earth’s orien- 
tation in inertial space, is March 21, 2000. In order to place the spacecraft at an initial 
terrestrial longitude of 75.07 deg (one of the stable longitudes), a true anomaly 0 of 
253.89 deg is used. 

These elements correspond to an initial position and velocity of 

f = -11698.237 7 — 40508.869 J + OK km 
v= 2.9547- 0.853 J + OK km/s 

The orbit control problem of geosynchronous satellites is a topic of continuing prac- 
tical interest. Detailed technical descriptions of standard north-south and east-west 
stationkeeping control techniques as well as more advanced orbit control concepts can 
be found in Refs. [11]-[13] and [18]-[20]. 

In the next section, we develop an attitude dynamics model of sun-pointing spacecraft 
in geosynchronous orbit for attitude control systems architecture design. 
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2.3 Rigid-Body Attitude Equations of Motion 

Consider a rigid body in a circular orbit. A local vertical and local horizontal (LVLH) 
reference frame A with its origin at the center of mass of an orbiting spacecraft has a 
set of unit vectors {ai, a 2 , a 3 } with d\ along the orbit direction, a 2 perpendicular to the 
orbit plane, and a 3 toward the earth, as illustrated in Figure 2.11. The angular velocity 
of A with respect to an inertial or Newtonian reference frame N is 

ujA/n = —na2 (131) 

where n is the constant orbital rate. The angular velocity of the body-fixed reference 
frame B with basis vectors {&i, f> 2 , & 3 } is then given by 

c o B/N = Co B/A + u A/N = Co B/A - na 2 (132) 

where uj b ! a is the angular velocity of B relative to A. 



Figure 2.11: A large space solar power satellite in geosynchronous orbit. 

The orientation of the body-fixed reference frame B with respect to the LVLH ref- 
erence frame A is in general described by the direction cosine matrix C as follows: 


b i 

b 2 



‘ C n 
C 2 i 

c 12 

c 22 

c 13 ' 
c 23 


" Si " 

(133) 

b 3 


C 31 

c 32 

C33 


. ^3 
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ai Cn C21 C31 b\ 

a 2 = C 12 C 22 C 32 b 2 (134) 

03 J L Cl 3 C 23 C 33 J 63 

When earth is modeled as a sphere with uniform mass distribution, the gravitational 
force acting on a small mass element with a mass of dm is given by 


df = ~ 


pRdm n(R c + p)dm 
~W~~ |£c + p1 3 


where p is the gravitational parameter of the earth, R and p are the position vectors 
of the small mass element from the earth’s center and the spacecraft’s mass center, 
respectively, and R c is the position vector of the spacecraft’s mass center from the 
earth’s center. 

The gravity-gradient torque about the spacecraft’s mass center is then expressed as: 


M p x df = —p, 


py_Rc 
l-Rc T p| 3 


and we have the following approximation 


I Rc + P\ — R c » 1 T 


2 (R c -p) , P 2 


Rr ~ 3 < 1 - 


3 (R c ■ P) 


higher-order terms 


where R c — \R C \ and p = |p|. By definition of the mass center, / pdm = 0; therefore, the 
gravity-gradient torque neglecting the higher-order terms can be written as 

M = % f(R c -p)(px R c )dm (138) 


This equation is further manipulated as follows: 


M = R c x J p{p- R c )dm 
= -j^Rc X J ppdm ■ R c 
= -%R C x [ [ p 2 Idm — J\- R c 


= -j^Rc x J P 2 1 dm ■ R c + j^R c x J ■ R c 
= RcxJ-Rc 
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since J = f(p 2 I — pp)dm and R c x I ■ R c = R c x R c = 0. The inertia dyadic of the 
spacecraft with respect to its mass center is denoted by J, and / represents the unit 
dyadic. 

Finally, the gravity-gradient torque is expressed in vector-dyadic form [17], [23] as: 

M = 3n 2 a 3 x J ■ a 3 (139) 


where n = \J /j / R\ is the orbital rate and a 3 = —R c /R c . 

In addition to the contribution to gravitational force (see Sec. 2.2.4), earth’s oblate- 
ness makes a contribution to gravity-gradient torque, shown in Ref. [23] to have a coeffi- 
cient of SplJ 2 RI/R 5 c . By comparing this to the coefficient above, 3/i/ R%, it is seen that at 
geosynchronous orbit the contribution of J 2 to gravity-gradient torque is approximately 
5 orders of magnitude less than the main term. 

The rotational equation of motion of a rigid body with an angular momentum H = 
J ■ uj b ! n in a circular orbit is then given by 


dH_ 

dt 




+ 0 B/N x H 

B 


M 


where {d/dt}^ indicates differentiation with respect to time in reference frame N, and 
{d/dt}s indicates differentiation with respect to time in reference frame B. The rela- 
tionship can be rewritten as 


J ■ £J + £) x J • 0 = 3n 2 a 3 x J ■ dz (140) 

where £0 = £J B ^ N , and note that £3 = {d£j/dt}N = {d£)/dt}B- 

Since £j, a 3 , and J can be expressed in terms of basis vectors of the body-fixed 
reference frame B as 

ui = uJ\bi + (*>2^2 T (141) 

az = Cizbi + C/ 23^2 + Czzbz (142) 

i=l j = 1 

the nonlinear equations of motion in matrix form become 
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To describe the orientation of the body-fixed reference frame B with respect to the 
LVLH reference frame A in terms of three Euler angles 0* (i = 1,2,3), consider the 
sequence of Ci(0i) <— C 3 (0 3 ) <— C 2 (0 2 ) from the LVLH reference frame A to a body- 
fixed reference frame B. For this rotational sequence, we have 


b i 


C n C 12 C \ 3 


’ ai “ 

b 2 

= 

C 2 1 c 22 c 23 



b s 


. c 3 1 c 32 c 33 


_ as 


c0 2 c0 3 s0 3 

— C 0 1 C 0 2 S 0 3 +S 0 l S 0 2 C 0 x C 0 3 
S 0 i C 0 2 S 0 3 + C 0 \ S 0 2 —S 0 iC 0 3 


— s 0 2 c 0 3 

C0i$0 2 S0 3 + S0iC0 2 
— S0\S0 2 S0 3 +C0iC0 2 



a\ 


a>2 


as 


where c 0, ■■ cos 0 t and s0i = sin 0*. 

And also we have the following kinematic differential equations: 


01 

0 2 

03 J 


1 


COS 0;>, 


COS 0 3 
0 
0 


— cos 0i sin 0 3 sin 0\ sin 0 3 
cos 0\ — sin 0i 

sin 0i cos 0 3 cos 0\ cos 0 3 



’ CJi ' 


' o ‘ 


&2 

+ 

n 


_ Cb?3 


0 


(143) 


The dynamical equations of motion about the body-fixed principal axes become 

JitO\ — ( J 2 — J 3 )ui 2 u) 3 = —3n 2 (J 2 — J 3 )C 23 C 33 (144) 

J 2 io 2 — (J 3 — J\)u 3 ui\ = -3n 2 (J 3 — Ji)C 33 C 13 (145) 

Js<jj 3 — (Ji — J 2 )u)\uj 2 = —3n 2 (Ji — J 2 )Ci 3 C 23 (146) 


where 


C\ 3 = — sin 0 2 cos 0 3 

C 23 = cos 0\ sin 0 2 sin 0 3 + sin 0i cos 0 2 

C 33 = — sin 0i sin 0 2 sin 0 3 + cos 0\ cos 0 2 

for the sequence of Ci(#i) <— C 3 (0 3 ) <— C 2 (0 2 ). 

One may linearize Eqs. (143)-(146) “about” an LVLH orientation while admitting a 
large pitch angle as follows. Assume 0\ and 0 3 remain small, allow 0 2 to be large, assume 
oj\ and uj 3 are small, and w 2 is equal to the sum of a small quantity and —n. Equations 
that are linear in the small quantities are 

Ji0i + (1 + 3cos 2 0 2 )n 2 (J 2 - J 3 )0 1 - n( Ji - J 2 + J 3 )0 3 
+ 3 (J 2 - J 3 )n 2 (sin 0 2 cos 0 2 )0 3 = u\ + di 
J 2 0 2 + 3n 2 ( Ji - J 3 ) sin 0 2 cos 0 2 = u 2 + d 2 
J ^03 T (1 T 3sin 2 ^ 2 )w 2 ( J 2 Ji)0 3 T n(</ J 2 ~\~ Js)0i 
+ 3 (J 2 - Ji)n 2 (sm0 2 cos0 2 )0i = u 3 + d 3 
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where u* and di are control and disturbance torques, respectively. 

For a quasi-inertially stabilized, sun-pointing SSP satellite in geosynchronous orbit 
with small body rates, a;* (i = 1, 2, 3), and small roll/yaw angles, 0i and 0 3 , the kinematic 
differential equations, (143), can be linearized in the small quantities: 

01 ~ uq 

02 OJ2 ~\~ Tt 

03 ~ c 03 

Finally, the equations of motion of a sun-pointing spacecraft with small roll and yaw 
angles can be found as 

J\9\ = — 3n 2 (02 — J3 )(cos 2 #2)#i — 3(J2 — J3)n 2 (sin02 cos 02)03 + U\ + di (147a) 
J 2 O 2 = — 3n 2 ( Ji — J 3 ) sin 0 2 cos O 2 + U 2 + d 2 (147b) 

■J 3 O 3 = - 3 n 2 (J 2 - Ji)(sin 2 0 2 )03 - 3(J 2 - Ji)n 2 (sin0 2 cos0 2 )0i + n 3 + d 3 (147c) 

The pitch angle relative to LVLH, 62, is not restricted to be small, but it may be 
regarded as a sum, 0 2 = nt + 86 2, where 40 2 is a small pitch attitude error. Kinematical 
and dynamical differential equations can then be made linear in the small quantities uq, 
u> 2 , u> 3 , #1, 402, and #3. For such a case, Eqs. (143) become 

0i « uq 
402 ~ Uq 
03 ~ (*>3 


and Eqs. (147) become 

Ji0i = — 3n 2 (J 2 — J 3 )[(cos 2 nt)0i + (sin nt cos nf)0 3 ] + Mi + d\ (148a) 

J2 89 2 = — 3n 2 ( Ji - J 3 )[(cos 2 nt — sin 2 nt)40 2 + sin nt cos nt] + u 2 + d 2 (148b) 

0303 = — 3 n 2 ( J 2 — Ji)[(sin 2 nt)0 3 + (sin nt cos nt)0i] + w 3 + d 3 (148c) 

where ( 01,402,03 ) are the small roll, pitch, and yaw attitude errors of a sun-pointing 

spacecraft, respectively. 

Equations (147) or (148) are the attitude equations of motion of the Abacus satellite 
for control design in the presence of the external disturbances, di, in units of N-m, 
modeled as: 


di ~ 12,000 — ll,900cosnf 

d 2 ~ 1, 200 

d 3 « —11, 900 sin nt 

However, ±20% uncertainties in this disturbance model as well as the inertia properties 
should be considered in control design. 
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2.4 Abacus Satellite Structural Models 


In Refs. [24]-[26], dynamics and control problems of large flexible platforms in space, 
such as the square Abacus platform, have been investigated. The flexible structure 
dynamics and control problem is a topic of continuing practical as well as theoretical 
interest. However, a significant control-structure interaction problem, possible for such 
very large Abacus platform (3.2 km by 3.2 km) with the lowest structural mode frequency 
of about 0.002 Hz, is avoided simply by designing an attitude control system with very 
low bandwidth (< orbit frequency of 1 x 10 -5 Hz). The proposed low-bandwidth attitude 
control system, however, utilizes a concept of cyclic-disturbance accommodation control 
to provide ±5 arcmin pointing of the Abacus platform in the presence of large external 
disturbances and dynamic modeling uncertainties. Consequently, the flexible structure 
control problem is not further investigated in this study, while a structural dynamic 
interaction problem with thermal distortion needs to be investigated in a future study. 

Various structural concepts for providing the required stiffness and rigidity of the 
Abacus platform are illustrated in Figure 2.12. Finite-element modeling of the baseline 
Abacus platform is illustrated in Figure 2.13 and the first three vibration modes are 
shown in Figure 2.14. Selected node locations for control analysis and design are shown 
in Figure 2.15. Typical pole-zero patterns of reduced-order transfer functions can be 
seen in Figure 2.16. Computer simulation results of a reduced-order structural model 
with the lowest 16 modes, confirm that the control-structure interaction problem can 
be simply avoided by the low-bandwidth attitude control system. Detailed technical 
discussions of the dynamics and control problems of flexible spacecraft can be found in 
the literature (e.g., see Refs. [17] and [27]), and thus the structural control problem of 
the Abacus satellite is not elaborated in this report. 
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ABACUS SUPPORT STRUCTURE CONFIGURATIONS 



Planar Array Prismatic Support Truss Tetrahedral 

(Built-Up Truss Beams) (modified configuration) Support Truss 



Figure 2.12: Abacus structural platform concepts (Courtesy of Tim Collins at NASA 
LaRC). 
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(16 arrays) 


MODIFIED ABACUS CONFIGURATION 
FINITE ELEMENT MODEL 


^ 3200m (80 arrays) ^ 



Square arrangement helps eliminate “weak” stiffness direction. 


Figure 2.13: Baseline Abacus finite element model (Courtesy of Tim Collins at NASA 
LaRC). 
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Planar Configuration, Thin Wall Struts 

(Free-Free Vibration Modes) 



Model, Mode 2, Mode 3, 

.0018 Hz .0026 Hz .0037 Hz 


2 r 



Figure 2.14: Baseline Abacus vibration modes (Courtesy of Tim Collins at NASA LaRC). 
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(16 arrays) 


Node Number Locations for Normal Modes Results 

(Nine Nodes Shown in Red) 

3200m (80 arrays) ^ 



Front View (Abaqus support truss in back) 


Figure 2.15: Selected FEM node locations for control analysis and design (Courtesy of 
Tim Collins at NASA LaRC). 
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Chapter 3 

Development of Abacus Control 
Systems Architecture 

3.1 Introduction to Control Systems Design 

This section provides a summary of the basic definitions and fundamentals in control 
systems design. It also provides the necessary background material for developing a 
control systems architecture for the Abacus satellite. Further detailed discussions of 
classical and modern control theory as applied to spacecraft control systems design can 
be found in Wie [17]. 

3.1.1 Feedback Control Systems 

Block diagram representations of a feedback control system are shown in Figure 3.1. 
Figure 3.1(a) is called a functional block diagram representation. Any physical system 
to be controlled is often called a plant. A set of differential or difference equations used to 
describe a physical system is called a mathematical model of the system. In the analysis 
and design of a feedback control system, we often deal with a mathematical model of 
the plant, not with the actual physical plant. Consequently, special care must be taken 
about uncertainties in the mathematical model because no mathematical model of a 
physical system is exact. 

A closed-loop feedback control system maintains a specified relationship between the 
actual output and the desired output (or the reference input) by using the difference of 
these, called the error signal. A control system in which the output has no effect on the 
control decision is called an open-loop control system. In a feedback control system, a 
controller, also called compensator or control logic, is to be designed to manipulate or 
process the error signal in order that certain specifications be satisfied in the presence 
of plant disturbances and sensor noise. In the analysis of control systems, we analyze 
the dynamical behavior or characteristics of the system under consideration. In the 
design or synthesis, we are concerned with designing a feedback control system so that 
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Disturbance 




Figure 3.1: Block diagram representations of a feedback control system. 


it achieves the desired system characteristics. 

A feedback control system can also be represented as in Figure 3.1(b) using transfer 
functions. In this figure, for simplicity, the actuator and sensor dynamics are neglected, 
and r{t) denotes the reference input, y(t) the plant output, G(s) the plant transfer 
function, K(s) the compensator, u(t) the control input, e(t) the error signal, w(t) the 
disturbance, d(t) the output disturbance, and n{t) a sensor noise. 

The output of this closed-loop system, neglecting the sensor noise n(t), can then be 
represented as 


V(s) 


r ( s ) + 3f) w f 8 

1 + K(s)G(s) { ’ 1 + K(s)G(s) { 


l + K{s)G{s) d ^ S 


( 1 ) 


where y(s ) = C[y(t)], r(s) = C[r(t)], w(s ) = £[«;(£)], and d(s) = £[d(t)\. In particular, 
the closed-loop transfer functions from d(s) and r(s) to y(s) are 


y(s) 

1 

= S( S ) 

( 2 ) 

d(s) 

" 1 + K(s)G(s) 

y(s) 

K(s)G(s) 

— T(s) 

( 3 ) 

r(s) 

" 1 + K(s)G(s) 
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and S(s) and Tfs) are called the sensitivity function and the complementary sensitivity 
function, respectively. Furthermore, we have the following relationship: 

S(s) + T(s) = 1 (4) 

The closed-loop characteristic equation is defined as 

1 + K(s)G(s) = 0 (5) 

and K(s)G(s) is called the loop transfer function. It is also called the open-loop transfer 
function in the literature. The importance of the loop transfer function cannot be 
overemphasized because it is used extensively in the analysis and design of closed-loop 
systems. The roots of the closed-loop characteristic equation are called the closed-loop 
poles. 

The error signal, ignoring the sensor noise n(t), is defined as 

e(t) = rft ) - yft) (6) 

and the steady-state error can be found as 

e ss = lim eft) = limse(s) (7) 

t— KX> S— ►() 


where e(s) = C[e[t)\, provided that eft) has a final value. For the system shown in 
Figure 3.1, ignoring w(s) and d(s), we have 



and 


Thus, it is required that 


sr(s) 

lim — 

3-0 l + K(s)G(s) 


lim K(s)G(s) 

S—> 0 


( 9 ) 

( 10 ) 


to have zero steady-state tracking error for a constant reference input command. 

A feedback control system is often characterized by its system type. The system type 
is defined as the number of poles of the loop transfer function K(s)G(s) at the origin. 
Therefore, a type 1 system has zero steady-state error for a constant reference input, a 
type 2 system has zero steady-state error for a constant or ramp reference input, and so 
forth. 


In order to reduce the effects of the disturbance, the magnitude of the loop transfer 
function K(s)G(s) must be large over the frequency band of the disturbance dft). For 
good command following at any frequency, the steady-state or D.C. gain must be large. 
In general, a fast transient response, good tracking accuracy, good disturbance rejection, 
and good sensitivity require a high loop gain over a wide band of frequencies. Because 
the high loop gain may degrade the overall system stability margins, proper tradeoffs 
between performance and stability are always necessary in practical control designs. 
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3.1.2 Classical Frequency-Domain Methods 
Root Locus Method 

One of classical control analysis and design techniques is the root locus method developed 
by Evans in 1950. In Evans’ root locus method, the closed-loop characteristic equation 
is described by 

l + KG(s) = 0 (11) 

where KG(s) denotes the loop transfer function, G(s) includes both the compensator 
transfer function and the plant transfer function, and K is called the overall loop gain. 
Note that the roots of the closed-loop characteristic equation are called the closed-loop 
poles. 

In Evans’ root locus plot, the poles and zeros of the loop transfer function KG(s ) 
are shown, where the poles are represented as cross, x, and zeros as circle, o. A root 
locus is then simply a plot of the closed-loop poles as the overall loop gain K is usually 
varied from 0 to oo. 

Using a root locus plot, one can easily determine a gain margin which is one of the 
most important measures of the relative stability of a feedback control system. A gain 
margin indicates how much the loop gain K can be increased or decreased from its chosen 
nominal value until the closed-loop system becomes unstable. For example, if the loop 
gain K can be increased by a factor of 2 until a root locus crosses the imaginary axis 
toward the right-half s plane, then the gain margin becomes 20 log 2 « +6 dB. In some 
cases of an open-loop unstable system, the closed-loop system may become unstable if 
the loop gain is decreased from its chosen nominal value. For example, if the gain can 
be decreased by a factor of 0.707 until the closed-loop system becomes unstable, then 
the (negative) gain margin is 20 log 0.707 ~ — 3 dB. The root locus method also allows 
the designer to properly select at least some of the closed-loop pole locations and thus 
control the transient response characteristics. 


Frequency-Response Methods 


Frequency-response analysis and synthesis methods are among the most commonly used 
techniques for feedback control system analysis and design, and they are based on the 
concept of frequency-response function. 

The frequency-response function is defined by the transfer function evaluated at s = 
ju)\ that is, the frequency response function of a transfer function G(s) is given by 

G(s) \ s=jlJJ = G(ju) = B&[G(ju))\ + j hn[G(ju))] = \G{juj)\e j4,{lJj) (12) 


where \G(joo)\ and <j>(u>) denote, respectively, the magnitude and phase of G(juj) defined 
as 


\G(ju)\ = ^/{Re[G(jo;)]} 2 + (Im[G(ju;)]} 2 
-i Im[G(ju;)] 


oj) = tan 


R e[G(ju)\ 
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For a given value of oj, G(joj) is a complex number. Thus, the frequency-response 
function G(ju) is a complex function of oj. Mathematically, the frequency response 
function is a mapping from the s plane to the G(joj) plane. The upper half of the 
ju- axis, which is a straight line, is mapped into the complex plane via mapping G(ju). 

One common method of displaying the frequency-response function is a polar plot 
(also called a Nyquist plot ) where the magnitude and phase angle of G{juj), or its real 
and imaginary parts, are plotted in a plane as the frequency to is varied. Another 
form of displaying G(joj) is to plot the magnitude of G(joj ) versus to and to plot the 
phase angle of G(juj) versus oj. In a Bode plot, the magnitude and phase angle are 
plotted with frequency on a logarithmic scale. Also, we often plot the magnitude of the 
frequency-response function in decibels (dB); that is, we plot 20 log \G(ju)\. A plot of 
the logarithmic magnitude in dB versus the phase angle for a frequency range of interest 
is called a Nichols plot. 

For a feedback control system shown in Figure 3.1, the loop transfer function, 
K(s)G(s) evaluated at s joj, is used extensively in the analysis and design of the 
system using frequency-response methods. The closed-loop frequency response func- 
tions defined as 


_ a/ ■ \ _ 

d(ju) J0J 1 + K(joj)G(joj) 
y(jw) = T( ■ \ = K(ju})G(ju}) 
r(ju) 1 + K(joj)G(joj) 


(13) 

(14) 


are also used in classical frequency-domain control systems design. 

One of the most important measures of the relative stability of a feedback control 
system is the gain and phase margins as defined as follows. 


Gain Margin. Given the loop transfer function K(s)G(s) of a feedback control 
system, the gain margin is defined to be the reciprocal of the magnitude \K(joj)G(joj)\ 
at the phase-crossover frequency at which the phase angle <f>(oj) is —180 deg; that is, the 
gain margin, denoted by g m , is defined as 


" \K(ju c )G(ju c )\ 


or 

g m = -20 log \K(joj c )G(joj c )\ dB (16) 

where oj c is the phase-crossover frequency. For a stable minimum-phase system, the 
gain margin indicates how much the gain can be increased before the closed-loop system 
becomes unstable. 


Phase Margin. The phase margin is the amount of additional phase lag at the 
gain-crossover frequency oj c at which \K(joj c )G(joj c )\ = 1 required to make the system 
unstable; that is, 

4>m = (t>[K{joj c )G{joj c )\ + 180° (17) 
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Although the gain and phase margins may be obtained directly from a Nyquist plot, 
they can also be determined from a Bode plot or a Nichols plot of the loop transfer 
function K{joj)G{joj). 

3.1.3 Classical PID Control Design 

The PID (proportional-integral-derivative) control logic is commonly used in most feed- 
back controllers. To illustrate the basic concept of the PID control, consider a cart of 
mass m on a frictionless horizontal surface, as shown in Figure 3.2(a). This so-called 
double integrator plant is described by 

my(t ) = u(t) + w(t) (18) 

where y is the output displacement of the cart, u is the input force acting on the cart, 
and w is a disturbance force. This system with a rigid-body mode is unstable, thus the 
system needs to be stabilized and the desired output is assumed to be zero. 

Assuming that the position and velocity of the system can be directly measured, 
consider a direct velocity and position feedback control logic expressed as: 

u(t) = -ky(t) - cy(t) (19) 


or 


u = —(k + cs)y 


where k and c are controller gains to be determined. The closed-loop system illustrated 
by Figure 3.2(b) is then described by 


my(t) + cy(t) + ky(t) = w(t ) 


which is, in fact, a mathematical representation of a mass-spring-damper system forced 
by an external disturbance w(t), as illustrated in Figure Figure 3.2(c). 

The closed-loop characteristic equation of the system shown in Figure 3.2 is 

ms 2 + cs + k = 0 


The control design task is to tune the “active damper” and “active spring” to meet 
given performance/stability specifications of the closed-loop system. Let co n and ( be 
the desired natural frequency and damping ratio of the closed-loop poles. Then the 
desired closed-loop characteristic equation becomes 

s 2 T 2 (^co n s T u>^ = 0 

and the controller gains c and k can be determined as 

c = 2m(u n 
k = 
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(20a) 

(20b) 




(a) Open-loop system 



(b) Closed-loop system with position and 
velocity feedback 



(c) Equivalent closed-loop system 
representation 

Figure 3.2: Control of a double integrator plant by direct velocity and position feedback. 
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The damping ratio £ is often selected as: 0.5 < £ < 0.707, and the natural frequency 
uj n is then considered as the bandwidth of the PD controller of a system with a rigid- 
body mode. For a unit-step disturbance, this closed-loop system with the PD controller 
results in a nonzero steady-state output y(o o) = l/k. However, the steady-state output 
error 2 /( 00 ) can be made small by designing a high-bandwidth control system. 

Consider the control problem of a double integrator plant with measurement of po- 
sition only. A common method of stabilizing the double integrator plant with noisy 
position measurement is to employ a phase-lead compensator of the form: 


u(s) 


—K 


T\S T 1 
T 2 s+ 1 


y(s) 


as illustrated in Figure 3.3(a). An equivalent closed-loop system can be represented 
using two springs and a damper as in Figure 3.3(b) and that 


K = h; 


Tt 


c{k\ T kf) 

hh ’ 


t 2 = 


c 

h 


For further details of designing a passive three-parameter isolator known as the D- 
Strut™ that can be modeled as Figure 3.3(b), see Davis, L. P., Cunningham, D., and 
Harrell, J., “Advanced 1.5 Hz Passive Viscous Isolation System,” Proceedings of AIAA 
Structures, Structural Dynamics, and Materials Conference, AIAA, Washington, DC, 
April 1994. 

In order to keep the cart at the desired position y = 0 at steady state in the presence 
of a constant disturbance, consider a PID controller of the form: 


u(t) = - K P y(t ) - AT/ J y(t)dt - K D y(t) (21) 

or 

u{s) = —[Kp + ~ + K d s\ y(s) 

In practical analog circuit implementation of a PID controller when y is not directly 
measured, differentiation is always preceded by a lowpass filter to reduce noise effects. 
It can be shown that for a constant disturbance, the closed-loop system with the PID 
controller, in fact, results in a zero steady-state output y( 00 ) = 0. 

The closed-loop characteristic equation can be found as 


ms T Kps T KpS T AT/ — 0 


and let the desired closed-loop characteristic equation be expressed as 

(s^ T 2£c u n s T u%)(s T 1/T 1 ) = 0 

where u n and £ denote, respectively, the natural frequency and damping ratio of the 
complex poles associated with the rigid-body mode and T is the time constant of the 
real pole associated with integral control. 
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(a) Closed-loop system with a phase-lead compensator 



(b) Equivalent closed-loop system representation using 
springs and a damper 


Figure 3.3: Control of a double integrator plant using a phase-lead compensator. 
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The PID controller gains can then be determined as 


(22a) 

(22b) 

(22c) 


i v- / 2 i > 

K P = m(u> n + -jr-) 

Ki=m^ 

1 

Kp = m(2(L0 n + — ) 


The time constant T of integral control is often selected as 




10 

C^n 


3.1.4 Digital PID Controller 

Consider a continuous-time PID controller represented as 


u(t) = —K P y(t) ~ Ki f y(t)dt - K D y(t ) 
Using Euler’s approximation of differentiation: 


1 — z 1 2 — 1 

T = T 2 


(23) 


we obtain an equivalent digital PID controller represented in z-domain transfer function 
form as: 

u = — Ik p + + K d — — — | y (24) 

This digital PID control logic can be implemented in a computer as follows: 


u(k) = -K P y(k) - Kju(k) - K d y(fe) ^ 


(25) 


where 

u(k) = u(k — 1) + Ty(k ) 

A single-axis block diagram representation of a digital control system of the Hubble 
Space Telescope is shown in Figure 3.4. As can be seen in this figure, the baseline digital 
control system of the Hubble Space Telescope, with a sampling period T = 0.025 sec 
and a computational delay of Td = 0.008 sec, is in fact a digital PID controller with a 
finite impulse response (FIR) filter in the rate loop. 
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HST DYNAMICS 


SOLAR ARRAY 



Figure 3.4: Simplified block diagram of the pitch-axis pointing control system of the 
Hubble Space Telescope [17], [28]. 
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3.1.5 Classical Gain-Phase Stabilization 


In the preceding sections, we have introduced the fundamentals of classical control. In 
this section, we present a classical gain-phase stabilization approach to compensator 
design, in particular, for a flexible spacecraft which has a rigid-body mode and lightly 
damped, oscillatory flexible modes. The approach allows the control designer to properly 
gain-phase stabilize each mode, one-by-one, resulting in a meaningful control design with 
physical insight. The classical gain-phase stabilization method is primarily restricted to 
the single-input single-output control problems, however. 

The classical concepts of gain-phase stabilization of a rigid-body and flexible modes 
can be summarized briefly as follows: 

1) Gain stabilization of a flexible mode provides attenuation of the control loop gain 
at the desired frequency to ensure stability regardless of the control loop phase 
uncertainty. A lightly damped, flexible mode is said to be gain stabilized if it is 
closed-loop stable for the selected loop gain, but it becomes unstable if the loop 
gain is raised or its passive damping reduced. Hence, a gain stabilized mode has a 
finite gain margin, but is closed-loop stable regardless of the phase uncertainty. 

2) Phase stabilization of a flexible mode provides the proper phase characteristics at 
the desired frequency to obtain a closed-loop damping that is greater than the 
passive damping of the mode. A lightly damped, flexible mode is said to be phase 
stabilized if it is closed-loop stable for arbitrarily small passive damping. Hence, a 
phase stabilized mode has a finite phase margin, but is closed-loop stable regardless 
of the loop gain uncertainty. 

3) A rigid-body or flexible mode is said to be gain-phase stabilized if it is closed-loop 
stable with finite gain and phase margins. 

When an actuator and a sensor are “colocated” on flexible structures in space, the 
rigid-body mode and all the flexible modes are said to be “stably interacting” with each 
other. For such a colocated case, position feedback with a phase-lead compensator or 
direct rate and position feedback can be used to stabilize all the flexible and rigid-body 
modes. Because all the modes are phase stabilized in this case, special care must be 
taken about the phase uncertainty from the control loop time delay and actuator /sensor 
dynamics. As frequency increases, the phase lag due to a time delay will eventually 
exceed the maximum phase lead of 90 degrees from the direct rate feedback. Thus, roll- 
off filtering (i.e., gain stabilization) of high-frequency modes is often needed to attenuate 
the control loop gain at frequencies above the control bandwidth. The selection of roll- 
off filter corner frequency depends on many factors. When a colocated actuator /sensor 
pair is used, the corner frequency is often selected between the primary flexible modes 
and the secondary flexible modes. An attempt to gain stabilize all the flexible modes 
should be avoided, unless the spacecraft or structures are nearly rigid. In practice, the 
actual phase uncertainty of the control loop must be taken into account for the proper 
tradeoff between phase stabilization and gain stabilization. 



When an actuator and a sensor are not colocated, the rigid body mode and some 
of the flexible modes are said to be “unstably interacting” with each other. Unless 
gain stabilization of all the flexible modes is possible for a low-bandwidth control, a 
proper combination of gain-phase stabilization is unavoidable. Gain stabilization of 
an unstably interacting flexible mode can be achieved only if that mode has a certain 
amount of passive damping. The larger the passive damping at a particular mode, the 
more conveniently it can be gain stabilized. Usually, gain stabilization is applied in 
order to stabilize high-frequency modes that have no significant effects on the overall 
performance. In practice, a structure has always a certain amount of passive damping, 
which allows for the convenient gain stabilization of such flexible modes. 

Notch filtering is a conventional way of suppressing an unwanted oscillatory signal 
in the control loop, resulting in gain stabilization of a particular flexible mode. The 
use of notch filtering ensures that the specific mode is not destabilized by feedback 
control; however, it does not introduce any active damping, which often results in too 
much “ringing” that may not be acceptable in certain cases. In general, roll-off of 
the control loop gain at frequencies above the control bandwidth is always needed to 
avoid destabilizing unmodeled high-frequency modes and to attenuate high-frequency 
noise, and it is often simply achieved by using a double-pole lowpass filter. To sharply 
attenuate a signal at high frequencies while affecting the magnitude and phase of the 
signal at low frequencies as little as possible, various high-order lowpass filters, such as 
Bessel, Butterworth, Chebyshev, or elliptical filters, are also used in feedback control 
systems, but mostly in open-loop signal processing. The common characteristic of these 
conventional filters is that they are minimum-phase filters. 

Although the last several decades have brought major developments in advanced 
control theory, the most usual approach to the design of practical control systems has 
been repetitive, trial-and-error synthesis using the root locus method by Evans and/or 
the frequency-domain methods by Bode, Nyquist, and Nichols. Classical control designs 
employ primarily a PID-type controller with notch and/or roll-off filtering. However, 
such classical control designs for a certain class of dynamical systems become difficult, 
especially, if a high control bandwidth is required in the presence of many closely spaced, 
unstably interacting, lightly damped modes with a wide range of parameter variations. 

For such case, the concept of generalized second-order filtering can be employed. The 
concept is a natural extension of the classical notch and phase lead/lag filtering, and it 
is based on various pole-zero patterns that can be realized from a second-order filter of 
the form 

s 2 /o^ + 2£ z s/w z + 1 
s 2 /ul + 2( p s/lu p + 1 

where u z , ( z , u p , and ( p are filter parameters to be properly selected. 

For different choices of the coefficients of this second-order filter, several well-known 
filters such as notch, bandpass, lowpass, highpass, phase-lead, and phase-lag filters can 
be realized. In addition to these minimum-phase filters, various nonminimum-phase 
filters can also be realized from this second-order filter [17]. 
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3.1.6 Persistent Disturbance Rejection 

A classical approach to disturbance accommodating control of dynamical systems in 
the presence of persistent or quasi-periodic disturbances is presented here. The method 
exploits the so-called internal model principle for asymptotic disturbance rejection. The 
concept of a disturbance rejection dipole is introduced from a classical control viewpoint. 

After successful stabilization of the rigid-body mode as well as any other unsta- 
bly interacting flexible modes, active disturbance rejection is then simply achieved by 
introducing into the feedback loop a model of the disturbance. A block diagram repre- 
sentation of a persistent disturbance rejection control system is shown in Figure 3.5. 

It is assumed that a persistent (or quasi-periodic) disturbance is represented as 


n 

w(t) = ^2 Ai sin(27r fit + 4>i) 

i— 1 


with unknown magnitudes and phases o x but known frequencies /*. Note that if, for 
example, f\ = 2/2 = • • • = nf n , then w(t) becomes a periodic disturbance. 

In general, the disturbance w(t) can be described by a Laplace transformation 


w(s) 


NUs) 

D w (s ) 


where N w (s ) is arbitrary as long as w(s) remains proper. The roots of D w (s) correspond 
to the frequencies at which the persistent excitation takes place. The inclusion of the 
disturbance model 1/D W inside the control loop is often referred to as the internal mod- 
eling of the disturbance. In classical design, the internal disturbance model is regarded 
as being part of the compensator as shown in Figure 3.5. The presence of 1/D W in the 
control loop results in the effective cancellation of the poles of w(s), provided that no 
root of D w (s) is a zero of the plant transfer function. This is shown in the following 
closed-loop transfer function: 


y(s) 


1 /D(s) 


-Ms) 


1 + N c (s)N(s)/D c (s)D w (s)D(s) 

D c (s)D m (s) Nwjs) 

D w (s)D c (s)D(s) + N c (s)N(s) D w (s ) 


(27) 


where we can see the cancellation of D w (s). 

The compensator can be viewed as a series of individual first-order or second-order 
filters as follows: 


N&) n NM 

D„(s) M DM 


Each filter is designed to perform a specific task, like the stabilization of a particular 
mode. In the same manner, a disturbance rejection filter can be designed that has a 
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Figure 3.5: Persistent disturbance rejection control system (transfer function descrip- 
tion). 
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proper transfer function and uses the internal disturbance model 1/D W . Thus a proper 
numerator is chosen in the compensator to go with the disturbance model as shown in 
Figure 3.5. The numerator is chosen to be of the same order as D w so that there is a 
zero for each pole of the disturbance model 1/D W . 

Although the asymptotic disturbance rejection based on the internal model principle 
has been well known, an interesting interpretation of the concept from a classical control 
viewpoint is presented here. Each pole-zero combination of the disturbance rejection 
filter 


n 


s 2 /co Zi 2 + 2 C Zi s/u Zi + 1 

s 2 AV + 1 


can be called a dipole, where C, Zi is included for generality. The filter thus consists of 
as many dipoles as there are frequency components in the persistent disturbance. The 
separation between the zero and the pole is generally referred to as the strength of the 
dipole. The strength of the dipole affects the settling time of the closed-loop system; 
in general, the larger the separation between the pole and zero of the filter the shorter 
the settling time is. This is caused by the position of the closed-loop eigenvalue corre- 
sponding to the filter dipole. As the strength of the dipole is increased, this eigenvalue 
is pushed farther to the left, speeding up the response time of the disturbance rejection. 
However, this separation influences the gain-phase characteristics of the system, because 
the dipole causes a certain amount of gain-phase changes in its neighborhood. More- 
over, at frequencies higher than the dipole there is a net gain increase or reduction. The 
magnitude of this gain increases with the separation between pole and zero. Therefore, 
as the strength of the dipole is changed to meet a chosen settling time the compensation 
must be readjusted. A compromise has to be reached often between the settling time 
and the stability of the compensated system. 

The internal model principle for persistent disturbance rejection is now incorporated 
with the standard state-space control design problem. Active disturbance rejection for 
the measured output y is to be achieved by introducing a model of the disturbance inside 
the control loop, therefore using again the concept of internal modeling, as illustrated 
in Figure 3.6. 

For example, consider a scalar disturbance d(t) with one or more frequency compo- 
nents represented as 

d(t) ■■ sin(o;*t + fa) 


with unknown magnitudes A* and phases fa but known frequencies uy. The disturbance 
rejection filter is then described by 


± d = A d x d + B d y 


(28) 
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Figure 3.6: Persistent disturbance rejection control system (state-space description). 


where x d is the state vector introduced by the disturbance model and, for example, 

0 10 0 ' 

-u\ 0 0 0 

0 0 0 1 ; ° d 

0 0 0 _ 

for a scalar output y(t) with d(t) of two frequency components. The disturbance rejection 
filter can include as many frequency components as the given disturbance, and is driven 
by the measured output y of the plant. This procedure is equivalent to the one used in 
the classical approach with the disturbance model now consisting of a state-space model. 
We now consider a plant described by the state-space equation: 

x p = A p x p + B p u + GpW (29a) 

y = C pX p + v + d (29b) 

where x p denotes the plant’s state vector, u the control input vector, w the process 
noise, v the measurement noise, and d the output equivalent persistent disturbance. 
Both w and v are assumed to be white noise processes with 

£^[w(t)w r (r)] = W 5(t — t ) 

-F[v(t)v r (r)] = V5(t — t) 

where W and V are the corresponding spectral density matrices. 
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In general, a compensator designed for this plant will consist of a regulator and an 
estimator which will approximate the states x p with estimated states x p using the infor- 
mation from the measured output y. The estimator which attempts to asymptotically 
reduce the error term e = x p — x p is given by 

x p = A p x p + B p u + L(y - C p x p ) 

= (A p - LCp)xp + BpU + Ly (30) 

where the term (y — C p Xp) represents the error between the output of the plant and the 
estimated output and L is the estimator gain matrix to be determined. 

The disturbance filter model described by Eq. (28) is then augmented to a plant 
described by Eq. (29) as follows: 


x = Ax + Bu + Gw 
y = Cx + v + d 

where 


X 

Xp 

*d _ 


A = 

Ap 

BdCp 

0 

Ad 

; b 

Bp 

0 

C 

Cp 

0 ] 

; G 

Gp 

0 





An estimated state feedback controller is then given as 


(31a) 

(31b) 


u = -Kx 

where x = \ x J xj J and the gain matrix K = [ K p K,j 1 is to be determined for 
the augmented system described by Eq. (31). 

As shown in Figure 3.6, however, x<j can be directly fed back as: 


u 


Kp K d 


Xp 


(32) 


since Xd is directly available from Eq. (28). 

An active disturbance rejection controller in state-space form is then given by 


£ 1 


_ x d J 



A p BpKp LCp 

0 


u = — 


Kp 



— BpK d 

Ad 



(33) 

(34) 


And the closed-loop system with w = v = d = 0is described as 


’ Xp " 


i 

-BpKp 

— BpKd 


" Xp " 

Xp 

= 

LCp 

Ap BpKp LCp 

— BpKd 


Xp 

. *d _ 


B d Cp 

0 

Ad 


_ *d _ 
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which can be modified using the error term, e = x p — x p , resulting in a partially decoupled 
system of equations, as follows: 


" x p “ 


A p BpKp 

— BjjK^ BpKj, 


Xp 

Xd 

= 

B d C p 

A d 0 


Xd 

e 


0 

0 A p — LC P _ 


e 


The closed-loop characteristic equation can then be written as 


si A p T BpKp 

—BjCp 

0 


BpKd 

si- A d 
0 


-BpKp 

0 

si A.p -t~ LCp 


0 


(35) 


The determinant in Eq. (35) is equal to the determinants of the diagonal submatrices 
multiplied together, one giving the regulator eigenvalues for the augmented system in- 
cluding the internal model, and the other giving the estimator eigenvalues for only the 
plant. Hence, we have shown that the separation principle for regulator and estimator 
holds for a closed-loop system even with an internal model for asymptotic disturbance 
rejection. 


3.1.7 Classical versus Modern Control Issues 

State-space approaches to control design are currently emphasized in the literature and 
more widely explored than classical methods. This arises from the convenience of obtain- 
ing a compensator for the whole system given one set of design parameters (e.g., given 
weighting matrices or desired closed-loop eigenvalues). In classical design, on the other 
hand, a compensator must be constructed piece by piece, or mode by mode. However, 
both classical and state-space methods have their drawbacks as well as advantages. All 
these methods require, nevertheless, a certain amount of trial and error. 

The question remains how to choose these parameters and what choice provides 
the “best” optimal design. The designer must find an acceptable set of parameters for a 
“good” optimal design. The use of state-space methods for control design usually results 
in a compensator of the same order as the system to be controlled. This means that 
for systems having several flexible modes, the compensator adds compensation even to 
modes that are stable and need no compensation. This may result in a complicated 
compensator design. 

The classical design is particularly convenient for the control of dynamical systems 
with well-separated modes. The concept of nonminimum-phase compensation also pro- 
vides an extremely convenient way of stabilizing unstably interacting flexible modes. 
The resulting compensator is usually of less order than the system to be controlled be- 
cause not all flexible modes in a structure tend to be destabilized by a reduced-order 
controller. A helpful characteristic of most flexible space structures is their inherent 
passive damping. This gives the designer the opportunity of phase stabilizing significant 
modes and to gain stabilize all other higher frequency modes which have less influence 
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Figure 3.7: Persistent-disturbance rejection control system for the ISS. 

on the structure. On the other hand, successive-mode-stabilization presents problems 
of its own, and a re-tuning of the compensated system becomes necessary. It is also 
noticed that reducing the damping in a frequency shaping filter reduces its influence on 
neighboring frequencies and it also reduces the phase lag at lower frequencies. However, 
reducing the damping of the filters increases the sensitivity of the phase stabilized modes 
to plant parameter uncertainties. 

Active disturbance rejection can be achieved in both the classical methods and state- 
space methods, with the introduction of an internal model of the disturbance into the 
feedback loop. The concept of internal modeling of the disturbance works as well with a 
classical transfer function description as with a state-space description. In the classical 
design, the internal modeling of the disturbance leads to the introduction of a disturbance 
rejection dipole, or filter, for each frequency component of the disturbance. In the state- 
space design the introduction of the internal model results in the addition of two states 
for each frequency component of the disturbance. 

Such a concept of persistent-disturbance rejection control has been successfully ap- 
plied to the International Space Station, as illustrated in Figure 3.7. Detailed control 
designs using a modern state-space control technique for the ISS, the Hubble Space 
Telescope, and large flexible structures can be found in [28]-[33]. 
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Table 3.1: Electric propulsion systems for the 1.2- GW Abacus satellite 


Thrust, T 

> 1 N 

Specific impulse, I sp = Tj {frig) 

> 5,000 sec 

Exhaust velocity, V e = I sp g 

> 49 km/s 

Total efficiency, 77 = P 0 /P* 

> 80% 

Power /thrust ratio, P/T 

< 30 kW/N 

Mass/power ratio 

< 5 kg/kW 

Total peak thrust 

200 N 

Total peak power 

6 MW 

Total average thrust 

80 N 

Total average power 

2.4 MW 

Number of 1-N thrusters 

> 500 

Total dry mass 

> 75,000 kg 

Propellant consumption 

85,000 kg/year 


Note: T = mV e , P 0 = |mE e 2 = \TV e , P 0 /T = jV e = ideal power/thrust ratio, Pi/T = 
%jV e , I sp = T/(fng ) = V e /g, V e = I sp g where g = 9.8 m/s 2 , fn is the exhaust mass flow 
rate, Pi is the input power, and P 0 is the output power. 

3.2 Control Systems Architecture 

The area-to-mass ratio of the Abacus satellite, A/m = 0.4 m 2 /kg, relatively large when 
compared to 0.02 m 2 /kg of typical geosynchronous communications satellites, is a key 
parameter characterizing the very large size of the Abacus satellite. If left uncontrolled, 
this can cause a cyclic drift in the longitude of the Abacus satellite of 2 deg, east and 
west. Thus, in addition to standard north-south and east- west stationkeeping maneu- 
vers for ±0.1 deg orbit position control, active control of the orbit eccentricity using 
electric thrusters with high specific impulse, I sp , becomes mandatory. Furthermore, 
continuous sun tracking of the Abacus satellite requires large control torques to counter 
various disturbance torques. A control systems architecture developed in this study uti- 
lizes properly distributed electric thrusters to counter, simultaneously, the cyclic pitch 
gravity-gradient torque and solar radiation pressure. 

Electric Propulsion Systems 

Basic characteristics of electric propulsion systems for the Abacus satellite are summa- 
rized in Table 3.1. 

Approximately 85,000 kg of propellant per year is required for simultaneous orbit, 
attitude, and structural control using 500 1-N electric propulsion thrusters with I sp = 
5,000 sec. The yearly propellant requirement is reduced to 21,000 kg if an I sp of 20,000 
sec can be achieved. As I sp is increased, the propellant mass decreases but the electric 
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Figure 3.8: A schematic illustration of the NSTAR 2.3-kW, 30-cm diameter ion thruster 
on Deep Space 1 Spacecraft (92-mN maximum thrust, specific impulse ranging from 
1,900 to 3,200 sec, 25 kW/N, overall efficiency of 45-65%). 


power requirement increases; consequently, the mass of solar arrays and power processing 
units increases. Based on 500 1-N thrusters, a mass/power ratio of 5 kg/kW, and a 
power/thrust ratio of 30 kW/N, the total dry mass (power processing units, thrusters, 
tanks, feed systems, etc.) of electric propulsion systems proposed for the Abacus satellite 
is estimated as 75,000 kg. 

A schematic illustration of the 2.3-kW, 30-cm diameter ion engine on the Deep Space 
1 spacecraft is given in Figure 3.8, which is formally known as NSTAR, for NASA Solar 
electric propulsion Technology Application Readiness system. The maximum thrust 
level is about 92 mN and throttling down is achieved by feeding less electricity and 
xenon propellant into the propulsion system. Specific impulse ranges from 1,900 sec at 
the minimum throttle level to 3,200 sec. 

In principle, an electric propulsion system employs electrical energy to accelerate 
ionized particles to extremely high velocities, giving a large total impulse for a small 
consumption of propellant. In contrast to standard propulsion, in which the products of 
chemical combustion are expelled from a rocket engine, ion propulsion is accomplished 
by giving a gas, such as xenon (which is like neon or helium, but heavier), an electrical 



charge and electrically accelerating the ionized gas to a speed of about 30 km/s. When 
xenon ions are emitted at such high speed as exhaust from a spacecraft, they push the 
spacecraft in the opposite direction. However, the exhaust gas from an ion thruster 
consists of large numbers of positive and negative ions that form an essentially neutral 
plasma beam extending for large distances in space. It seems that little is known yet 
about the long-term effect of such an extensive plasma on geosynchronous satellites. 

Orbit, Attitude, and Structural Control System 

A control systems architecture developed in this study is shown in Figure 3.9. The 
proposed control systems utilize properly distributed ion thrusters to counter, simul- 
taneously, the cyclic pitch gravity-gradient torque, the secular roll torque caused by 
cm-cp offset and solar pressure, the cyclic roll/yaw microwave radiation torque, and the 
solar pressure force whose average value is 60 N. A control-structure interaction prob- 
lem of the Abacus platform with the lowest structural mode frequency of 0.002 Hz is 
avoided simply by designing an attitude control system with very low bandwidth (< 
orbit frequency). However, the proposed low-bandwidth attitude control system utilizes 
a concept of cyclic-disturbance accommodating control to provide ±5 arcmin pointing 
of the Abacus platform in the presence of large external disturbances and dynamic mod- 
eling uncertainties. High-bandwidth, colocated direct velocity feedback, active dampers 
may need to be properly distributed over the platform. 

Placement of approximately 500 1-N electric propulsion thrusters at 12 different 
locations is illustrated in Figure 3.10. In contrast to a typical placement of thrusters 
at the four corners, e.g., employed for the 1979 SSPS reference system, the proposed 
placement shown in Figure 3.10 minimi zes roll/pitch thruster couplings as well as the 
excitation of platform out-of-plane bending modes. A minim u m of 500 ion engines of 1- 
N thrust level are required for simultaneous attitude and stationkeeping control. When 
reliability, lifetime, duty cycle, lower thrust level, and redundancy of ion engines are 
considered, this number will increase significantly. 


3.3 Control System Simulation Results 

Computer simulation results of a case with initial attitude errors of 10 deg in the presence 
various dynamic modeling uncertainties (e.g., ±20 % uncertainties in moments and 
products of inertia, center-of-mass location, and principal axes, etc.), but without cyclic- 
disturbance rejection control, are shown in Figures 3.11-3.15. It can be seen that the 
pointing performance is not acceptable. 

Control simulation results of a case with 10-deg initial attitude errors in the presence 
various dynamic modeling uncertainties (e.g., ±20 % uncertainties for inertia, cm loca- 
tion, and principal axes, etc.), and with additional cyclic-disturbance rejection control, 
are shown in Figures 3.16-3.20. The proposed low-bandwidth attitude control system 
that utilizes the concept of cyclic-disturbance accommodation control satisfies the ±5 
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arcmin pointing requirement of the Abacus platform in the presence of large external 
disturbances and dynamic modeling uncertainties. Proper roll/pitch thruster firings 
needed for simultaneous eccentricity and roll/pitch attitude control can be seen in Fig- 
ure 3.19. Nearly linear control forces are generated by on-off modulation of individual 
1-N thrusters, as can be seen in this figure. The total thrusting force from the roll/pitch 
thrusters #1 through #4 nearly counters the 60-N solar pressure force. 
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Feedforward Control 


Torque Command Gravity -Gradient T orque 

u 2c= 3n 2 (J 1 -J 3 )(sin2nt)/2 -3n 2 (Jj- J 3 ) (sin 2 0 2 )/2 



Figure 3.9: An integrated orbit, attitude, and structural control system architecture 
employing electric propulsion thrusters. 
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Thrust force direction 



Roll: 1/3 Pitch: 2/4 Yaw: 5/6/7/8 
Orbit Eccentricity, Roll/Pitch Control: 1/3, 2/4 

E/W and Yaw Control: 9/10/11/12 
N/S and Yaw Control : 5/6/7/8 


Figure 3.10: Placement of a minimum of 500 1-N electric propulsion thrusters at 12 
different locations, with 100 thrusters each at locations #2 and #4. (Note: In contrast 
to a typical placement of thrusters at the four corners, e.g., employed for the 1979 
SSPS reference system, the proposed placement of roll/pitch thrusters at locations #1 
through #4 minimizes roll/pitch thruster couplings as well as the excitation of platform 
out-of-plane bending modes.) 
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Figure 3.11: Simulation results without cyclic-disturbance rejection control. 


103 







Pitch Rate (deg/s) 


x 10' 4 



Time (Orbits) 

Figure 3.12: Simulation results without cyclic-disturbance rejection control (continued). 
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Figure 3.13: Simulation results without cyclic-disturbance rejection control (continued). 
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Figure 3.14: Simulation results without cyclic-disturbance rejection control (continued). 
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Figure 3.15: Simulation results without cyclic-disturbance rejection control (continued). 
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Figure 3.16: Simulation results with cyclic-disturbance rejection control. 
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Figure 3.17: Simulation results with cyclic-disturbance rejection control (continued). 
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Figure 3.18: Simulation results with cyclic-disturbance rejection control (continued). 
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Figure 3.19: Simulation results with cyclic-disturbance rejection control (continued). 
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Figure 3.20: Simulation results with cyclic- disturbance rejection control (continued). 
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Chapter 4 

Conclusions and Recommendations 


4.1 Summary of Study Results 

The major objective of this study was to develop advanced concepts for controlling or- 
bit, attitude, and structural motions of very large Space Solar Power Satellites (SSPS) 
in geosynchronous orbit. This study focused on the 1.2-GW “Abacus” SSPS concept 
characterized by a square (3.2 x 3.2 km) solar array platform, a 500-m diameter mi- 
crowave beam transmitting antenna, and an earth-tracking reflector (500 x 700 m). For 
this baseline Abacus SSPS configuration, we derived and analyzed a complete set of 
mathematical models, including external disturbances such as solar radiation pressure, 
microwave radiation, gravity-gradient torque, and other orbit perturbation effects. An 
integrated orbit, attitude, and structural control systems architecture developed for the 
Abacus satellite employs properly distributed, 500 1-N electric propulsion thrusters. 

Despite the importance of the cyclic pitch gravity-gradient torque, this study shows 
that the solar pressure force is considerably more detrimental to control of the Abacus 
satellite (and other large SSPS) because of an area-to-mass ratio that is very large 
compared to contemporary, higher-density spacecraft. 

A key parameter that characterizes the sensitivity of a satellite to solar radiation 
pressure is the area-to-mass ratio, A/m ; the value of A/m for the Abacus satellite is 0.4 
m 2 /kg, which is relatively large when compared to 0.02 m 2 /kg for typical geosynchronous 
communications satellites. Solar radiation pressure causes a cyclic drift in the longitude 
of the Abacus satellite of 2 deg, east and west. Consequently, in addition to standard 
north/south and east /west stationkeeping maneuvers for ±0.1 deg orbit position control, 
active control of the orbit eccentricity using electric thrusters becomes nearly mandatory. 
Furthermore, continuous sun tracking of the Abacus platform requires large control 
torques to counter various disturbance torques. 

The proposed control system architecture utilizes properly distributed ion thrusters 
to counter, simultaneously, the cyclic pitch gravity-gradient torque, the secular roll 
torque caused by center of mass - center of pressure offset and solar pressure, the cyclic 
roll/yaw microwave radiation torque, and the solar pressure force whose average value 
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is 60 N. In contrast to a typical placement of thrusters at the four corners, e.g., em- 
ployed for the 1979 SSPS reference system, the proposed placement shown in Figure 
3.10 minimizes roll/pitch thruster couplings as well as the excitation of platform out-of- 
plane bending modes. A control-structure interaction problem of the Abacus platform 
with the lowest structural mode frequency of 0.002 Hz is avoided simply by designing 
an attitude control system with very low bandwidth (< orbit frequency). However, the 
proposed low-bandwidth attitude control system utilizes a concept of cyclic disturbance 
accommodation to provide ±5 arcmin pointing of the Abacus platform in the presence of 
large external disturbances and dynamic modeling uncertainties. Approximately 85,000 
kg of propellant per year is required for simultaneous orbit, attitude, and structural 
control using 500 1-N electric propulsion thrusters with a specific impulse of 5000 sec. 
Only 21,000 kg of propellant per year is required if electric propulsion thrusters with a 
specific impulse of 20,000 sec can be developed. As I sp is increased, the propellant mass 
decreases but the electric power requirement increases; consequently, the mass of solar 
arrays and power processing units increases. 

The total dry mass (power processing units, thrusters, tanks, feed systems, etc.) of 
electric propulsion systems for the Abacus satellite is estimated as 75,000 kg based on 
a minimum of 500 1-N thrusters and a mass/power ratio of 5 kg/kW. The peak power 
requirement is estimated as 6 MW based on the total peak thrust requirement of 200 N 
and a power /thrust ratio of 30 kW/N. 


4.2 Recommendations for Future Research 

The baseline control system architecture developed for the Abacus satellite requires 
a m inimum of 500 ion engines of 1-N thrust level. The capability of present electric 
thrusters are orders of magnitude below that required for the Abacus satellite. If the 
xenon fueled, 1-kW level, off-the-shelf ion engines available today, are to be employed, 
the number of thrusters would be increased to 15,000. The actual total number of ion 
engines will further increase significantly when we consider the ion engine’s lifetime, relia- 
bility, duty cycle, redundancy, etc. Consequently, a 30-kW, 1-N level electric propulsion 
thruster with a specific impulse greater than 5,000 sec needs to be developed for the 
Abacus satellite if excessively large number of thrusters are to be avoided. 

Several high-power electric propulsion systems are currently under development. For 
example, the NASA T-220 10-kW Hall thruster recently completed a 1,000-hr life test. 
This high-power (over 5 kW) Hall thruster provides 500 mN of thrust at a specific 
impulse of 2,450 sec and 59% total efficiency. Dual-mode Hall thrusters, which can 
operate in either high-thrust mode or high-/ sp mode for efficient propellant usage, are 
also being developed. 

The exhaust gas from an electric propulsion system forms an essentially neutral 
plasma beam extending for large distances in space. Because little is known yet about 
the long-term effect of an extensive plasma on geosynchronous satellites with regard 
to communications, solar cell degradation, etc, the use of lightweight, space-assembled 
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Table 4.1: Technology advancement needs for the Abacus SSPS 



Current 

Enabling 

Electric Thrusters 

3 kW, 100 mN 

30 kW, 1 N 


I sp = 3, 000 sec 

I sp > 5, 000 sec 


(5,000-10,000 thrusters) 

(500-1,000 thrusters) 

CMGs 

20 N-m-s/kg 

2,000 N-m-s/kg 


5,000 N-m-s/unit 

500,000 N-m-s/unit 


(500,000 CMGs) 

(5,000 CMGs) 

Space- Assembled 


66,000 N-m-s/kg 

Momentum Wheels 


4 xlO 8 N-m-s/unit 

(300-m diameter) 


(5-10 MWs) 


large-diameter momentum wheels may also be considered as an option for the Abacus 
satellite; therefore, these devices warrant further study. The electric thrusters, CMGs, 
and momentum wheels are compared in Table 4.1 in terms of their technology advance- 
ment needs. It is emphasized that both electrical propulsion and momentum wheel 
technologies require significant advancement to support the development of large SSPS. 

Despite the huge size and low structural frequencies of the Abacus satellite, the 
control-structure interaction problem appears to be a tractable one because the tight 
pointing control requirement can be met even with a control bandwidth that is much 
lower than the lowest structural frequency. However, further detailed study needs to 
be performed for achieving the required 5-arcmin microwave beam pointing accuracy 
in the presence of transmitter /reflector-coupled structural dynamics, Abacus platform 
thermal distortion and vibrations, hardware constraints, and other short-term impulsive 
disturbances. 

Although the rotating reflector concept of the Abacus satellite eliminates massive 
rotary joint and slip rings of the 1979 SSPS reference concept, the transmitter fixed to 
the Abacus platform results in unnecessarily tight pointing requirements imposed on the 
platform. Further system-level tradeoffs will be required for the microwave-transmitting 
antenna design, such as whether or not to gimbal it with respect to the platform, use 
mechanical or electronic beam steering, and so forth. 

The following research topics of practical importance in the areas of dynamics and 
control of large flexible space platforms also need further detailed investigation to support 
the development of large SSPS. 

• Thermal distortion and structural vibrations due to solar heating 

• Structural distortion due to gravity-gradientloading 

• Simultaneous eccentricity and longitude control 

• Attitude control during the solar eclipses 
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• Orbit and attitude control during assembly 

• Attitude and orbit determination problem 

• Reflector tracking and pointing control problem 

• Microwave beam pointing analysis and simulation 

• Space-assembled, large-diameter momentum wheels 

• Electric propulsion systems for both orbit transfer and on-orbit control 

• Backup chemical propulsion systems for attitude and orbit control 
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Appendix A 

Simulation of Orbital Motion 


A.l Introduction 

Numerical simulations of orbital motion, the results of which are presented in Chapters 
1 and 2, employ the algorithms described briefly in what follows. 

Encke’s method, as described in Sec. 9.4 of Ref. [16], and in Sec. 9.3, of Ref. [2], lies 
at the heart of a MATLAB/SIMULINK computer program used to integrate dynamical 
and kinematical equations governing relative translational motion of two bodies. 

This appendix begins with a brief description of the general relationship for two-body 
motion, then provides an overview of Encke’s method and how it is carried out in the 
computer program, and ends with a presentation of the expressions used in computing 
the various contributions to the perturbing forces exerted on the two bodies. 


A. 2 Two-Body Motion 

As discussed in Chapter 2, the relative orbital motion of two bodies is described by 

?+ = f = fs-fp (1) 

where f is the position vector from the mass center P* of a planet P to the mass center 
B* of a body B, r is the magnitude of r, f indicates the second derivative of f with 
respect to time t in an inertial or Newtonian reference frame N, and /i = G(mp + rrip), 
where G is the universal gravitational constant, mp is the mass of P, and rrip is the 
mass of B. 

If P were a sphere with uniform mass distribution, or a particle, and if B were 
a particle, then the gravitational force exerted by P on B would be given by g = 
— Gmpmpf/r 3 . The force exerted by B on P would be simply — g. The vector fp 
represents the resultant force per unit mass acting on B, other than g/mp', fp represents 
the resultant force per unit mass acting on P, other than —g/mp. 
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When / is as large or larger than yur/r 3 , integration of Eq. (1) is advisable and is 
referred to as Cowell’s method. On the other hand, when / is small in comparison to 
/ir/r 3 Cowell’s method can be disadvantageous in terms of numerical efficiency, and a 
different strategy known as Encke’s method may be preferred. 


A. 3 Encke’s Method 


The method of Encke requires the solution of ordinary differential equations governing 
the behavior of 5, 

6 = f — p (2) 

where p represents the solution of Eq. (1) with / 0; the path traced out by p is a 

conic section, known as the osculating orbit. The orbit described by r is the actual or 
true orbit of B about P, which differs from the osculating orbit whenever / does not 
vanish. 

The behavior of 5 is governed by Eq. (9.27) of Ref. [16], 

1 (3) 


where S indicates the second derivative of S with respect to time t in N, and p is the 
magnitude of p. The function / of q is given by 


m = q 


3 + 3g + g 2 

1 + (1 + q)i 


(4) 


where q is defined as 

qthlzp. (5) 

The values of 5 and 5 are both zero at the beginning of each simulation, and also 
following orbit rectification, or the point at which the osculating position and veloc- 
ity, p and p, are made equal to the true position and velocity, f and r, respectively. 

Rectification is performed when, as suggested in Ref. [2], (5*5) > 0.01(p • p) ' . 

The osculating orbit is determined as a function of time using initial values for p and 
p (which change with each rectification), together with Battin’s universal formulae for 
conic orbits according to Eqs. (3.33) and (4.84), and the relationships given in Prob. 4- 
21 of Ref. [1]. Use of the universal formulae requires a generalized anomaly x, obtained 
by Newtonian iteration as set forth in Eq. (4.4-15) of Ref. [2], or at the top of p. 219 in 
Ref. [1]; iteration is terminated when the time associated with x through the generalized 
form of Kepler’s equation [Eq. (4.81), Ref. [1]], is within 1 x 10 -4 sec of the simulation 
time t. 
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Six scalar, first order, ordinary differential equations corresponding to the second 
order vector Eq. (3) are integrated using a variable step, Runge-Kutta 4-5 scheme, with 
relative and absolute error tolerances set to 1 x 10 -8 . The true position and velocity, f 
and r, are used to calculate classical orbital elements a, e, i, Q, lo, and M according to 
the material in Secs. 2.3 and 2.4 of Ref. [2], and Secs. 3.3 and 4.3 of Ref. [1]. 


A. 4 Contributions to the Perturbing Force 

In the case of geosynchronous satellites the perturbing force per unit mass / receives 
significant contributions from the gravitational attraction of the Sun and Moon, Earth’s 
tesseral gravitational harmonics of degree 2 and orders 0 and 2, and solar radiation 
pressure, as discussed in Sec. 2.2. The remainder of this section contains the expressions 
employed in the computer program for these contributions, denoted respectively as f s , 
fm, f 2,0, A, 2 , and fr, such that 

f— fs + fm + f2,0 + A, 2 + fr (6) 


A. 4.1 Solar and Lunar Gravitational Attraction 


The gravitational force per unit mass exerted by the Sun on P is given by /x s r s /r 8 , where 
!i s is the product of G and the Sun’s mass, f s is the position vector from P* to the Sun’s 
mass center, and r s is the magnitude of r s . Likewise, the gravitational force per unit 
mass exerted by the Sun on B is given by p, s (f s — f)/\f s — r| 3 . Therefore, 


Hs{r s - f) n a f 3 



(7) 


When f'is small in comparison to f s , numerical difficulties can be encountered in the 
evaluation of the right hand member of Eq. (7); therefore, an alternate form of f s is 
used, as suggested in Eq. (8.61) of Ref. [1]: 


fs 



[f + ] 


( 8 ) 


where 



f • (r — 2r s ) 
f s • f s 


(9) 


The position vector f s from Earth’s mass center (actually, the Earth-Moon barycen- 
ter) to the Sun’s mass center, projected onto geocentric-equatorial directions and referred 
to the ecliptic of date, is obtained as a function of t with the formulae and numerical 
values given on p. E4 of Ref. [3]. 

Similarly, the contribution of lunar gravitational attraction to / is given by 


Pm ( r) 

| fm - f \ 3 


f^mX n 


( 10 ) 
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where n m is the product of G and the Moon’s mass, f m is the position vector from P* 
to the Moon’s mass center, and r m is the magnitude of f m . Numerical difficulties are 
avoided by using the expression 

fm = -■ y^P— [f+ f(q m )f m ] (11) 

\> m I | 


where 


Qrn 


a r • (f - 2 f m ) 


m * n 


( 12 ) 


The position vector f m from Earth’s mass center to the Moon’s mass center, projected 
onto geocentric-equatorial directions, and referred to the mean equator and equinox of 
date, is obtained as a function of t with the algorithm set forth on p. D46 of Ref. [3]. 


A. 4. 2 Tesseral Harmonics 


The computer program makes use of Eq. (12) of Ref. [4] to account for the gravitational 
harmonics of P, for any degree n and order m; in the simulations performed for this 
study, n and m are limited to 2. Numerical values of the gravitational coefficients, 
gravitational parameter of Earth, and mean equatorial radius, are those of the Goddard 
Earth Model T1 as reported in Ref. [5]. 

Earth’s oblateness is represented by a zonal harmonic of degree 2 and order 0, and 
is responsible for precessions in a satellite’s orbit plane and argument of perigee. The 
contribution of this harmonic to the force per unit mass exerted by P on B is given in 
Eq. (45) of Ref. [4] (also Prob. 3.7b in Ref. [6]) as 


f _ T R© 2 

/ 2,0 - ~ 


3 sin (be , 3 + ^(1 — 5 sin 2 <f>)-\ 
2 r 


(13) 


where n® is the gravitational parameter of the Earth, the product of G and the Earth’s 
mass; R© is the mean equatorial radius of the Earth (6378.137 km), r is the magnitude 
of r, and <f> is the geocentric latitude of B. Unit vector ejj is fixed in the Earth in the 
direction of the north polar axis. 

The contribution of oblateness to the force per unit mass exerted by B on P is given 
by —tub /2,0/mp, and the contribution of oblateness to / is thus [1 + (mfi/mp)]/^. In 
the case of the SSP orbiting Earth = 25 x 10 6 kg and mp = 5.98 x 10 24 kg, so 
mp/mp = 4x 10 -18 , which can be neglected in comparison to 1; therefore, the entire 
contribution of oblateness to / is essentially equal to f 2 ,o- 

The contribution f 2> 1 of the tesseral harmonic of degree 2 and order 1 vanishes because 
the harmonic coefficients S 2 , 1 and 62,1 are both zero. The harmonic of degree 2 and order 
2 can cause the longitude of of a geosynchronous spacecraft to drift; from Eq. (12) of 
Ref. [4] the contribution to the force per unit mass exerted by P on B is given by 


A, 2 — 


£«©R© 2 ( C 2}2 C 2 + S 2>2 S 2 


42,363 — (sin </>4 2 ,3 + 5^2,2) — 
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(14) 


+ 2^2,2 [(^ 2 , 2^1 + S' 2 , 2 <Sl)ei + (52,2^1 — C 2t2 S\ )e 2 \ | 

where unit vectors e x and e 2 are fixed in the Earth: e\ lies in the equatorial plane 
parallel to a line intersecting Earth’s geometric center and the Greenwich meridian, and 
e 2 = e 3 x e x . 

Equations (6) and (7) of Ref. [4] indicate that the required derived Legendre poly- 
nomials are A 2 ,2 = 3 and A 2> 3 = 0. In addition, Eqs. (9) and (10) of Ref. [4] show 
that 

Si = f • e 2 = r cos ^ sin A, C\ = f • e\ = r cos ^ cos A (15) 

S 2 = 2(rcos</>) 2 sinAcosA, C 2 = (rcos0) 2 (cos 2 A — sin 2 A) (16) 

where A is the geographic longitude of B measured eastward from the Greenwich merid- 
ian. Therefore, 


f 2> 2 = | — 15 cos 2 ^>[62,2 (cos 2 A — sin 2 A) + 2S 2>2 sin A cos A] ^ 

-I- 6 cos <t> [((72,2 cos A + S 2>2 sin A)ei + (S 2i2 cos A - C 2>2 sin A)e 2 ] | (17) 

As in the case of f 2i o, ms/mp is neglected in comparison to 1, and f 2>2 thus constitutes 
the entire contribution of the present harmonic to /. 


A. 4. 3 Solar Radiation Pressure 

The force per unit mass of solar radiation pressure exerted on B is given by —C{f s — 
f)/(m,B\f s — f|) where C is a constant, 60 N. We neglect the solar radiation pressure 
exerted on the Earth, and write 


J = Cjf-fs) 
m B \f s ~f\ 


(18) 
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